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Wc present a detailed implementation of two parallel bispectrum estimation methods which can be 
applied to general non-separable primordial and CMB bispectra. The method exploits bispectrum 
mode decompositions on the tctraliodral domain of allowed wavenumber or multipole values, using 
both separable basis functions and related orthonormal modes. We provide concrete examples 
of such modes constructed from symmetrised tetrahcdral polynomials, demonstrating the rapid 
convergence of expansions describing nonseparable bispectra. We use these modes to create rapid 
and robust pipelines for generating simulated CMB maps of high resolution [l > 2000) given an 
arbitrary primordial power spectrum and bispectrum or an arbitrary late-time CMB angular power 
spectrum and bispectrum. By extracting coefficients for the same separable basis functions from 
an observational map, we are able to present an efficient /nl estimator for a given theoretical 
model with a nonseparable bispectrum. The estimator has two manifestations, comparing the 
theoretical and observed coefficients at either primordial or late times, thus encompassing a wider 
range of models, including secondary anisotropics and leasing as well as active models, such as 
cosmic strings. We provide examples and validation of both /nl estimation methods by direct 
comparison with simulations in a WMAP-realistic context. In addition, we demonstrate how the 
full primordial and CMB bispectrum can be extracted from observational maps using these mode 
expansions, irrespective of the theoretical model under study. We also propose a universal definition 
of the bispectrum parameter Fnl, so that the integrated bispectrum on the observational domain 
can be more consistently compared between theoretical models. We obtain WMAP5 estimates of 
/nl for the equilateral model from both our primordial and late-time estimators which are consistent 
with each other, as well as results already published in the literature. These general bispectrum 
estimation methods should prove useful for nonGaussianity analysis with the Planck satellite data, 
as well as in other contexts. 

I. INTRODUCTION 

Standard inflationary scenarios predict the Universe to be close to flat with primordial curvature pertur- 
bations which are nearly scale-invariant and Gaussian. All these predictions are in very good accord with 
cosmic microwave background (CMB) and large-scale structure measurements, such as those provided by 
WMAP and SDSS. Despite this remarkable agreement, present observations are not able to completely rule 
out alternatives to inflation, nor to effectively discriminate among the vast number of different inflationary 
models that have been proposed. However, almost all such quantitative comparisons derive from inferred 
measurements of the primordial two-point correlator or power spectrum P{k) from {(,(,), where C, is the 
curvature perturbation. If wc wish to subject inflation to more stringent tests and to distinguish between 
competing models then perhaps the best prospects are offered by studying nonGaussianity, that is, the 
higher order correlators beyond the power spectrum. The three-point correlator of the CMB or bispectrum 
Bij^i^ is a projection on the sky of the evolved primordial bispectrum B{ki, k2, k^) arising from (CCC)) con- 
sisting of contributions from triangle configurations with sidclcngths given by the wavenumbers ki,k2,k^. 
The bispectrum has attracted most attention in the literature to date and its study is usually simplified 
to the characterization of a single nonlinearity parameter /nl, which schematically is given by the ratio 
hi,^B{k,k,k)/P{kf. 

Standard inflation, that is, single field slow-roll inflation, predicts a very small bispectrum with /nl ~ 0.01 
[1, 2], possessing a characteristic scale-invariant local shape. (This local shape is dominated by squeezed 
triangle conflgurations, that is, those for which one side is much smaller than the others, e.g. ki <C /c2, k^.) In 
fact, such a low signal would be undetectable even by an ideal noiseless CMB experiment, because it is below 
the level of NG contamination expected from secondary anisotropics /nl ~ C'(l)- However, measurement 
of a significantly larger primordial /nl ^ 1 would have profound consequences because it would signal the 



2 




Figure 1: Flow chart for the two general estimator methodologies described and implemented in this article using complete 
separable mode expansions. Note the overall redundancy which assists estimator validation and the independence of the 
extraction of expansion coefficients from theory q:„ (cycle f ) and data /3„ (cycle 2). Explanations for the schematic equations 
can be found in the main text. 



need for new physics during inflation or even a paradigm shift away from it. Present measurements of this 
local /nl are equivocal with the WMAP team reporting [3] 

/nl = 51 ± 60 (95%) (1) 

and with other teams obtaining higher [4] (WMAP3) or equivalent values [5, 6], while with improved 
WMAP5 noise analysis a lower value was found /nl = 38 it 42, but at a similar 2a significance [7]. The 
Planck satellite experiment is expected to markedly improve precision measurements with A/nl = 5 or 
better [8]. 

Further motivation for the study of the bispectrum comes from the prospect of distinguishing alternative 
more complex models of inflation which can produce nonGaussianity with potentially observable amplitudes 
/nl ^ li but also in a variety of different bispectrum shapes, that is, with the nonGaussian signal peaked 
for different triangle configurations of wavevectors. To date only special separable bispectrum shapes 
have been constrained by CMB data, that is, those that can be expressed (schematically) in the form 
B{ki, k2, k^) = X (ki)Y (k2) Z (ks) , or else can be accurately approximated in this manner. All CMB analysis, 
such as those quoted above for the local shape (1), exploits this separability to reduce the dimensionality of 
the required integrations and summations to bring them to a tractable form. The separable approach reduces 
the problem from one of 0{l^^^) operations to a manageable 0{l^^^) [9]. Other examples of meaningful 
constraints on separable bispectrum shapes using WMAP5 data include those for the equilateral shape [3] 
and another shape 'orthogonal' to both equilateral and local [10]. Despite these three shapes being a good 
approximation to non-Gaussianity from a number of classes of inflation models, they are not exhaustive 
in their coverage of known primordial models [11], nor other types of late-time non-Gaussianity, such as 
that from cosmic strings [12, 13]; they cannot be expected to be, given the functional degrees of freedom 
available. Bringing observations to bear on this much broader class of cosmological models, therefore, is 
the primary motivation for this paper. 

In a previous paper [14], we described a general approach to the estimation of non-separable CMB bis- 
pectra. The method has developed out of the first direct calculations of the reduced CMB bispectrum bi^i^i^ 
which surveyed a wide variety of non-separable primordial models, revealing smooth coherent patterns of 
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acoustic peaks on the tetrahedral domain of allowed multipole values. Since the hj^i^ could be well repre- 
sented using a limited number of bins, we could similarly decompose them into orthogonal mode functions 
which converged in relatively short mode expansions [11]. Here, we describe the detailed implementation of 
these methods in a comprehensive dual approach to estimating bispectrum parameters which is illustrated 
in fig. 1. We present concrete examples of separable basis functions Q„ (symmetrised tetrahedral polynomi- 
als) and corresponding orthonormal modes TZn on the domain of allowed wavenumbers ^1,^2,^3; these are 
then deployed within a more general mode expansion methodology. In the first primordial implementation, 
we decompose an arbitrary non-separable shape S using separable basis functions with coefficients a„. This 
expansion can be used for a fast calculation of the full CMB bispectrum Bi^i^i^ (section III), as well as lead- 
ing to a robust method for generating simulated maps from a given power spectrum P{k) and bispectrum 
B{ki, k2, ks) (section IV). Our main emphasis here, however, is on a primordial estimator for /nl which is 
achieved by a confrontation between theory, represented by the a„ coefficients, and a set of observational 
coefficients found by extracting the same modes from the observed CMB map (section III). Examples of 
simulated maps and recovery of the input /nl are given in section V in a WMAP-realistic context. 

In the second and parallel late-time implementation (see fig. 1), we assume the theoretical CMB bispec- 
trum -B/iia/s is calculated already from the primordial shape [11] or because it is a late-time effect ranging 
from secondary anisotropies through to fluctuations induced by cosmic strings. A separable mode expansion 
of Bi^i^i^ allows for a simpler and more direct approach to /nl estimation, as well as simulated map genera- 
tion, in a wider variety of scenarios. Here, as well as primordial models we consider the antithetical example 
of cosmic strings. These two estimator methods are complementary with each having distinct advantages 
depending on the properties and generation mechanism of the non-Gaussianity under investigation. They 
provide independent validation in situations where both are applicable. 

It remains to point out recent and related developments, especially those by colleagues in Planck Working 
Group 4 (NonGaussianity) . To date most primordial shapes have been assumed to be scale-invariant, but 
in ref. [15] some deviations from the local shape were considered in developing a more general approach. 
Spherical Mexican wavelets, using a limited number of scales, were employed in ref. [5, 16] to estimate /nl 
for the local shape with WMAP3 data, providing a constraint consistent with (1). Similar work has been 
achieved for needlets with corresponding constraints [6] ) , again essentially tailoring the method to the local 
template using local shape map simulations. Another approach to a late-time estimator has also exploited 
the smoothness of the reduced CMB bispectrum by using a limited number of multipole bins [17]. The 
method was tested for the local shape using map simulations, and emphasised Planck forecasts investigating 
the pattern of acoustic peaks in the local model. We shall discuss here how these late time approaches - 
whether wavelets, bins or other alternatives - fall within the general mode expansion methodology outlined 
previously [14] and can be applied, in principle, to explore nonseparable primordial models beyond local 
nonGaussianity. We point out in the implementation presented here, however, that direct estimation of the 
CMB bispectrum can be achieved without reference to the calculated bispectrum for a particular model 
and without relying on corresponding CMB map simulations. 



In this section we will review some basic definitions and mathematical formulae that will be used through- 
out the rest of the paper. Our work will be concerned with the analysis of the three-point function induced 
by a NG primordial gravitational potential $(k) in the CMB temperature fluctuation field. Temperature 
anisotropies are represented using the aim coefficients of a spherical harmonic decomposition of the cosmic 
microwave sky, 



II. THE CMB BISPECTRUM AND /nl ESTIMATION 



A. Relation between primordial and CMB bispectra 
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The primordial potential $ is imprinted on the CMB mutipoles ai„i by a convolution with transfer functions 
A;(fc) representing the linear perturbation evolution, through the integral 

^AKA;)$(k)ll„(k). (2) 

The CMB bispectrum is the three point correlator of the aim, so substituting we obtain 

) (3) 

($(ki)$(k2)$(k3)) Yi^mAk,)Yi^m2{h)Yi^m,{ks) , (5) 

where fci = |ki|, /u2 = |k2| and /cs = |k3|. Here, we define the primordial bispectrum as 

($(ki)$(k2)$(k3)) = {27TfB^{ki, fc2, h) ,5(ki + k2 + ka) , (6) 

where the delta function enforces the triangle condition, that is, the constraint imposed by translational 
invariance that wavevectors in Fourier space must close to form a triangle, ki + k2 + ka = 0. We replace 
the delta function in (6) with its exponential integral form, substitute this into equation (3) and integrate 
out the angular parts of the three kj integrals in the usual manner to yield 

^™m3 = 0) J X^dx J dkidk2dk3{kik2h)^B^{ki,k2,k3)Ai^{h)Ai^{k2)Ai^{k3) 

X jhikix)ji^{k2x)ji^{k3x) J (i0^y,^^^(x)yi,^,(x)yi3„3(x). (7) 

The last integral over the angular part of x is known as the Gaunt integral which can be expressed in terms 
of Wigner-3j symbols as 

(2/i + l)(2/2 + l)(2/3 + l) (lil2h\( h h h \ 



47r Y y y mi m2 ma 

Given that most theories we shall consider are assumed to be isotropic, it is usual to work with the angle- 
averaged bispectrum, 

Bhhh = XI ( ri\ rn2 ma ) ^"'^"^^'^'^^'"^"^ams) • (9) 

or the even more convenient reduced bispectrum which removes the geometric factors associated with the 
Gaunt integral. 

The reduced bispectrum from (3) then takes the much simpler form 

Kkh = J x^dx J dkidk2dk3 {kik2k3f S$(A;i, fe, fca) 

X Ai^(ki) Ai^{k2) AiXk3) ji^ikix) 3i^{k2x) ji^{k3x) . (11) 

Here, it is important to note that the Gaunt integral in (10) encodes several constraints on the angle 
averaged bispectrum Bijj^ which are no longer transparent in the reduced bispectrum bijj^. These are. 
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first, that the sum of the three muhipoles must be even and, secondly, that the ^j's satisfy the triangle 
condition, analogously to the wavenumbers ki. For wavenumbers, the triangle condition is enforced through 
the x-integral over the three spherical Bessel functions jiihx) which evaluates to zero if the ki's cannot 
form a triangle, whereas in multipole space it is enforced by the angular integration dflx over the spherical 
harmonics Yj-m; in (8). Appreciating the origin of these constraints is important when we later consider the 
separability of the reduced bispectrum expression (7). 



B. Separable primordial shapes and CMB bispectrum solutions 

Given that the primordial power spectrum is very nearly scale-invariant, it is expected that the bispectrum 
will behave similarly. In order to bring the bispectrum to a scale-invariant form we have to appropriately 
eliminate a scaling which naturally arises in (6). This is usually achieved by multiplying through by the 
factor {kik2k3)^ appearing in (11) and defining a primordial shape function as 

S{ki, k2, ks) = ^{kik2k3fB^{ku k2, ks) , (12) 

where A'^ is a normalisation factor which is often taken such that for equal ki the shape function has unit 
value S{k, k, k) = 1. (This normalisation is also used for /nl, but it only strictly applies for scale-invariance 
and, in any case, leads to inconsistent comparisons between different models, as we shall discuss in section 
IV.) We thus characterise scale-invariant models in terms of an overall amplitude, parametrised by /nl, 
and their transverse shape, described by S{ki, k2, k^) on a triangular slice with ki + k2 + k^ = const. [18]. 
This leaves a two-dimensional space on which it is most elegant to use the two independent variables a, f3 
[14, 19] 

a = {k2 — kz) /k , (3 = {k — ki)/k, where k = ^{ki + k2 + k^) = const. , (13) 

with the following domains < /3 < 1 and — (1 — /3) < a < 1 — (3. For scale-dependent models with a 
non-trivial variation in k, the full three-dimensional dependence on the ki must be retained. In terms of 
the shape function (12), the reduced bispectrum (11) can be rewritten as 

Khh = ^(^) / ^^'^^ / dkidk2dk3S{ki,k2,k3)Ai^{ki)Ai^{k2)Ai^{k3)ji^{kix)ji^{k2x)ji^^^^^ 

The simplest possible shape function is the constant model 

Siki,k2,k3) = l, (15) 

for which a large-angle analytic solution for the reduced bispectrum was presented in ref. [11], 
Al 1 



LConst $ 

hhh - 27Ar(2/i + l)(2Z2 + l) (2/3 + 1) 



+ 



li + l2 + h + 3 h + h + l. 



(/ < 200) . (16) 



Here, we take the Sachs-Wolfe approximation that A; (k) = ((to — T^gc) k) for / <C 200 and exploit the 
manifest separability of the expression (14) to perform the one-dimensional ki integrations individually. 
The more general constant solution does not have an analytic solution for I > 200, for the reason that 
the transfer functions cannot be expressed in a simple form, but it can be evaluated numerically from the 
expression 

Cws = wj ^^dxli,i^)Ii,{x)Ii,{x), where Ii{x) = ^J dk Ai{k) ji{kx) . (17) 

The large-angle solution (16) is an important benchmark with which to compare the shape of late-time 
CMB bispectra from other models hi^i^i^ (note the scaling) and, additionally, it has some further recent 
physical motivation [20]. 
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The most studied scale-invariant shape function is the local model, 

where, in the second line, we also allow for power spectra which are nearly scale invariant, defined by 
($(k)$*(k')) = (27r)^P(A;)5(k-k') with P{k) ~ k~^. Using the Sachs- Wolfe approximation again, this has 
the corresponding large-angle analytic solutions 

.local ^ 2A| / 1 1 1 \ 

'^'^'3 277r2 + iMh + 1) hih + l)l3{k + 1) ^3(^3 + iMh + l)J ^ ' 

Here, we see that the divergences for the squeezed triangles {ki <C k2, k^...) in the primordial shape (18) are 
also reflected in ^^-king it a much less useful for relative comparison than the constant model (16). 

It is straightforward, in principle, to calculate the full bispectrum from the separable expressions arising 
ftom (18), 

bul = I ^'dx [ai,{x)0i^{x)Pi^{x) + (2 perms)] , (20) 

where the separated integrals analogous to (17) become 

ai{x) = ^Jdkk^ Ai{k) jiikx) , f3i{x) = ljdk k^P{k) Ai{k) ji{kx) . (21) 

We note that these highly oscillatory integrals must be evaluated numerically with considerable care. 
The separable equilateral shape has also received a great deal of attention with [18] 

, , X (fci + k2- k^){k2 + fc3 - fci)(fc3 + fci - fca) 
S{ki,k2,k3) = — — (22) 

fclfe2fc3 



= -2- 



^2 

y — h (2 perms) 



^2^3 



+ 



+ (5 perms) 

k2 



This is a much more regular shape than local (18) with the signal dominated by equilateral triangle con- 
figurations fci ~ ^2 ~ ^3 (the apparent divergence of the local shape in the second term cancels against the 
third). There is no simple large-angle analytic solution known for the equilateral model, unlike (17) and 
(20). In order to calculate the full equilateral bispectrum we evaluate the simplified expression 

Cms = / x'dx{2di,SiJi, + [oci.^iAs + (2 perms)] + [A.7i.5^3 + (5 perms)] } , (23) 

where cti, /3; are given in (21) and Ji, Si are defined by 

-Tiix) = ^Jdkk^ P{ky/^Ai{k) jiikx) , Mx) = lJ dk k^ P{kf/^Ai{k) jiikx) . (24) 

The equilateral shape is not derived directly from a physical model, but was chosen phenomenologically 
as a good separable approximation to specific models including the non-local part of Maldacena's original 
shape [1], as well as non-canonical cases such as higher derivative models [21] and DBI inflation [22] (for a 
review of single- field inflation shapes, see e.g. ref. [23]). These shapes are, in general, non-separable from the 
perspective of the integral (14). Here, we give a specific shape example for a model with higher derivative 
operators (which is also identical to DBI inflation): 

= k^k2k^ik^ + ^2 + ks)^ ( ^ ^ " ^^'''^^ ^ ^ ^^"^^'^^ " ^^^'^^^'^ ) ' ^^^^ 
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Not only is the equilateral shape (22) an excellent approximation to (25), a full Fisher matrix analysis of the 
respective CMB bispectra has shown they are 99% correlated out to /max < 2000 [11]. However, a simple 
separable approximation is not necessarily available for arbitary primordial shapes, nor is a particular 
separable representation necessarily convenient from a calculational perspective (as we shall discuss in 
section V for the equilateral case above). In ref. [11], we reviewed models currently proposed in the 
literature showing that families of CMB bispectra arising from non-separable shapes, such as feature and 
flattened models, are largely independent of the separable models currently constrained observationally (see 
also discussion of a 'cosine' shape correlator in ref. [18]). The independence of two shapes S and S' can be 
calculated from the integral [11] 

Fe{S,S')= [ S{ki,k2,k3)S'{ki,k2,h)uJe{kl,k2,k3)dVk, (26) 

where we choose the weight to be 

yj^ki,k2,ks) = -—^-—-, (27) 
ki + k2 + k3 

reflecting the scaling we see in the CMB correlator we meet in the next section. The shape correlator is 
then defined by 

^ ' ' ^yF{s,s)F{s',s') ■ ^ ^ 

By way of further illustration of the need to move beyond simple separable primordial shape functions, 
we present the late-time CMB bispectrum predicted analytically for cosmic strings [13] 



^string ^ 
I1I2I3 



2I3 50L I V 500 



ill -H- ll) {^ + A)J^ erf(0.3C/3) + 2 perms 



{I < 2000) , (29) 



where Imin = min(^i, ^2, ^3), h = min(500, Imin), C = min(l/500, 1/lmin) and 

L = c^i(qq + iiq + qq)-i(if + 4 + ii) . (so) 

Here, A ~ (SttG//)^ is a model dependent amplitude with Gfj, = /u/rripj measuring the string tension fj, 
relative to the Planck scale. The cutoffs around I ~ 500 in (29) are associated with the string correlation 
length at decoupling (perturbations with / > 500 can only be causally seeded after last scattering). (For 
the original small angle solution valid for I ^ 2000, see ref. [12, 13].) Here, the non-separable nature and 
very different scaling of the string CMB bispectrum are clear from a comparison with (19). Moreover, given 
the late-time origin of this signal from string metric perturbations, the modulating effect of acoustic peaks 
from the transfer functions is absent. 



C. Estimators for /nl and related correlators 



The main purpose of this non-Gaussian CMB analysis is to measure the CMB bispectrum induced by non- 
Gaussianities in the primordial gravitational potential, the link being given by equation (11). Unfortunately, 
the bispectrum signal is too weak to measure individual multipoles directly, so to compare theory with 
observation we must use an estimator which sums over the available multipoles. An estimator can be 

thought of as performing a least squares fit of the bispectrum predicted by theory {ai^mj^ai^ni-i^hms) to 
the bispectrum actually obtained from observations 

«^mx<m. ^ats- Ignoring sky cuts and inhomogeneous 
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noise, the estimator is weighted with the expected signal variance from Ci and written in the simple form 

(o/imi 0-hm2 C'i,r'«3 ) ""^mi '^ifms ^hma 



N ^ 



-I ^obs ^obs ^obs 

/ y ='mim2m3 ^hhh ■ V"-* / 

where we have used (10) and the Gaunt integral is given in (8) and N is the usual normalisation factor, 



We note from the second line of (31) that, for a given theoretical model, we need only calculate the reduced 
bispcctrum hi^i^i^ rather than the much more challenging full bispectrum, {ai^^^ai^m^ai^m^ . 

The above estimator has been shown to be optimal [24] for general bispectra in the limit where the 
non-Gaussianity is small and the observed map is free of instrument noise and foreground contamination. 
Of course, this is an idealised case and we need to consider taking into account the effect of sky cuts and 
inhomogeneous noise, which was considered in some detail in refs [25, 26]. In the more general case the 
optimal estimator takes the form: 



S = 



^E(ii\ X j^/iMsX (33) 



C-^a°^') (c-^a"^') (c-^a"^') + Cr^ (c-'^a"^' 



where the covariance matrix C is now non-diagonal due to mode-mode coupling introduced by the mask 
and anisotropic noise. Moreover, due to the breaking of isotropy, an additional term linear in the aim lias 
now to be added in order to maintain the optimality of the estimator [24]. In the ideal case one can easily 
see that the linear term is proportional to a monopole, while the covariance matrix is diagonal and equal 
to 1/C;, thus reproducing the initial formula (31). 
In this paper we will follow the approach of [27] and approximate the estimator (33) as 



I y ^m^^'^ms Khh f obs obs _ g™ \ obs 



where the tilde denotes modification to include experimental effects. The normalisation becomes 

^=fskyY.^f^^ (35) 

with the C/'s and h^^i^i^ now incorporating beam and noise effects through 

Ci = biCi + Ni and hiM3 = KKKhM3- (36) 

Here, hi is the beam transfer function, Ni the noise power spectrum, f^^y the fraction of the sky remaining 
after application of the mask and C^"^^^ is the covariance matrix calculated from Gaussian simulations. 
In what follows, it will be clear from the context whether beams, noise and masks are being incorporated 
in the analysis, so for simplicity we shall continue with the original estimator notation (31). 

The estimator (31) also naturally defines a correlator for testing whether two competing bispectra could 
be differentiated by an ideal experiment. Replacing the observed bispectrum with one calculated from a 
competing theory we have. 
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where now the normalisation N is defined as follows, 



N 



n/2 



3 



While this late time correlator is the best measme of whether two CMB bispectra are truly independent, 
it requires a full calculation of the CMB bispectrum which is time consuming in general. In [11] we 
determined that for the majority of models the shape correlator (28) introduced earlier is sufHcent to 
determine independence. 

An inspection of equations (31,33) shows that a brute force numerical implementation of the optimal 
estimator above would take 0(t^^^ operations. This means an implementation is not feasible for the 
angular resolutions achieved by present and forthcoming datasets (e.g. in the signal dominated regime we 

have Iraax < 500 for WMAP and I 

max ^ 2000 for Planck). However, as initially shown in ref. [9], if a specific 
theoretical bispectrum can be written in separable form as B{ki, k2, k^) = X (ki)Y {k2) Z (k^) then the 
computational cost of the algorithm can be reduced to 0{l'}^^^) operations, making the estimation tractable 
even at very high angular resolutions. This establishes the fact that separability is a crucial property for 
realistic data analysis, even though it is not generic for well-motivated inflationary and other models. As we 
have seen, the usual solution adopted has been to approximate the primordial non-separable shape under 
study using a separable form that is highly correlated with the original. This kind of approach requires a 
case-by-case analysis of all non-separable bispectra arising from different models and an educated "guess" of 
a good separable approximation, the close correlation of which must be verified numerically before moving 
on to the real analysis. Besides being impractical, this can also prove to be extremely difficult in specific 
cases. The aim of this work is then to find a completely general mathematical framework to "separate" 
shapes, both primordial and late-time, and thus build a general pipeline for /nl estimation and simulation 
of non-Gaussian CMB maps, that can be applied to any shape of interest. 



III. BISPECTRUM MODE DECOMPOSITION 



Our goal is to represent arbitrary non-separable primordial bispectra B{ki,k2,k?,) or CMB bispectra 
bij^i^ on their respective wavenumber or multipole domains using a rapidly convergent mode expansion 
[14]. Moreover, we need to achieve this in a separable manner, making tractable the three-dimensional 
integrals required for bispectrum estimation (14) by breaking them down into products of one-dimensional 
integrals. In particular, this means that we wish to expand an arbitrary non-separable primordial shape 
function as 

S{ki,k2, ^3) = XI E E Ipi^i) Qrih) Qsih) , (39) 

p r s 

where the qp are appropriate basis mode functions which are convergent and complete, that is, they span 
the space of all functions on the bispectrum wavenumber (or multipole) domain. In what follows below, 
we present one pathway for efficiently achieving this objective in stages. First, we create examples of one- 
dimensional mode functions qp{ki) in the /ci-direction which are orthogonal and well-behaved over the full 
wavenumber (or multipole) domain. We then construct three-dimensional products of these mode functions 
Qp{ki)qr{k2)Qp{k3) — *• Qn creating a complete basis for all possible bispectra on the given domain. Finally, 
by orthonormalising these product basis functions Qn TZn, we obtain a rapid and convenient method for 
calculating the relevant expansion coefficients aprs in (39). The subsequent discussion and implementation 
of general primordial and CMB bispectrum estimators, as well as map-making methods, is then built 
around these mode functions qp, Qn, and TZn- Here, we use bounded symmetric polynomials as a concrete 
and working implementation of this methodology, and we defer discussion about other possible basis mode 
functions which have been investigated to the end of the section. 
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Figure 2: Tetrahedral domain ('tetrapyd') for allowed multipole values / for the CMB bispectrum bi-^i^i^ or, with wavenumbers 
k for the primordial bispectrum B{ki, k2, fca))- The regular tetrahedral region defined up to the equilateral slice h + h + h < 
2imax = 2L (shaded brown) contains two thirds of the overall volume. The rest of the domain is given by the regular triangular 
pyramid on top which fills the volume to the corner of the encompassing cube defined by h, h, h < L. An origami tetrapyd is 
also shown (right) with folding instructions. 



In Fourier space, the primordial bispectrum B{ki, k2, k^) is defined when the three wavevectors ki, k2, ka 
close to form a triangle ki + k2 + ks = 0. Since each such triangle is uniquely defined by the lengths of its 
sides ki = |ki|, ^2 = |k2|, ks = jksl, we only require wavenumbers in the bispectrum argument. In terms of 
these three wavenumbers, the triangle condition restricts the allowed combinations into a tetrahedral region 
defined by 

ki < k2 + ks for ki > k2, k^, or k2 < ki + k^ for k2 > ki, k^, or k^ < ki + k2 for ^3 > ki, k2 ■ (40) 

This region forms a regular tetrahedron if we impose the restriction that ki + k2 + k^ < 2/cinax) however, 
it is more natural to extend the domain out to values given by a maximum wavenumber in each direction 
ki,k2,k3 < /cmax- This extension is motivated by issues both of separability and observation. The allowed 
domain Vr is then a hexahedron formed by the intersection of a tetrahedron and a cube. It can be obtained 
from a regular tetrahedron (two-thirds of the total volume) by gluing on top a regular triangular pyramid 
constructed from the corner of the cube (as illustrated in fig. 2). For brevity, let us denote this asymmetric 
triangular bipyramid as a tetrapyd, from the merger of a tetrahedron and a pyramid. Of course, bispectrum 
symmetries are such that it is only necessary to use one sixth of this domain, but aesthetics and intuition 
are helped by keeping the full domain while making a restriction to symmetrised functions. 

We will frequently need to integrate functions f{ki,k2,k^) over the tetrapyd domain (40), which for 
brevity we will denote as V-r with the integration given explicitly by 



A. Tetrahedral domain and weight functions 




(41) 



{ to'' H'' H-^ FWdzdx dy i;-! FWdzdy dx+ 

+li/2 FWdzdydx + J^-y /;_^ FWdzdxdy} 



where K = /Cmax) w{ki,k2,k^) is an appropriate weight function, and we have made the transformation 
X = ki/K, y = k2/K, x = k^/K with F{x,y,z) = f{Kx,Ky,Kz) and W{x,y,z) = w[Kx,Ky,Kz). For 
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Figure 3: The reduced CMB bispectra for the equilateral model (left) and the local model (right) plotted on the tetrahedral 
region shown in figure 42 (from [11]). Several density contours are illustrated (light blue positive and magenta negative) and 
biii2i3 is normalised by scaling relative to the constant Sachs- Wolfe solution (16) bY^°2h /b'l'^rSs- ^o^e the acoustic peaks induced 
by the transfer functions and the centre weighting for the equilateral model, contrasting with the corner-weighting for the local 
case [14]. 

integrals over the product of two functions / and g we can define their inner product (/, g) = T[fg], 
essentially defining a Hilbert space of possible shape functions in the domain (40). The total volume of the 
tetrapyd domain is given by T[l] = Initially, for the sake of simplicity, on the primordial wavenumber 

domain we will restrict attention to unit sidelength K = 1 and weight w = 1. 

We note that it is important to incorporate a weight function for a variety of reasons. For example, the 
primordial shape function S{ki, k2, k^) can be shown to possess a nearly linear scaling with respect to the 
CMB bispectrum estimator; on the multipole domain wij^i^ is non-trivial. A fairly close correspondence 
between the two can be obtained using w{ki,k2,k^) ~ 1/(^1 + k2 + k^) [11] which explains its choice in 
the shape correlator (28). The choice of weight function also affects mode expansion convergence and for 
certain shapes it may be convenient to eliminate dependencies by rescaling with a separable fuction. For 
the shapes we consider here however, this is not necessary. 

When analysing the CMB bispectrum it is particularly important to extend the tetrahedral domain to 
include multipoles in the top pyramidal region shown in fig. 2. In principle, this pyramid contains 33% of the 
triple /i^2^3 combinations available in the observational data, e.g. with Planck out to hjhjh ^ 2000. The 
tetrapyd domain for the reduced bispectrum bij^h becomes the discrete {hjh^h} combinations satisfying 

max ) 

^1 < ^2 + ^3 for li > I2, I3, + cyclic perms. , (42) 
^1 + ^2 + ^3 = 2n , n G N . 

In fig. 3 we illustrate contrasting bispectra on this domain for the equilateral and local models (here with 

/max = 2000). 

In multipole space, we will be primarily dealing with a summation over all possible {hjh, h} combinations 
in the estimator (31) or the closely related correlator (37). The appropriate weight function in the sum is 
then 

^IM. = h (2^1 + l)(2/2 + mh + 1) ( \ \ 0' , (43) 
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where we note that the third condition in (42) arises as a selection rule from the Wigner-3j symbol. 
Despite the discrete origin of the function wijj^, like the reduced bispectrum bijj^, it varies smoothly. 
It is particularly uniform on cross-sectional slices h + h + h = 2L, except for a finite rise very close 
to the boundaries. While the Wigner-3j symbols are easily calculable (especially in the rrij = case 
when performed in advance for a look-up table), it is more convenient to work in the continuum limit 
^Iil2l3~^ w{h, hi h) when considering domains with large /max- To achieve this we take the exact expression 
in terms of factorials (for even combinations with li + I2 + I3 = 21 , / G N), 



/i^2/3\ , {2l-2h)\{2l-2l2)\{2l-2k)\ l\ 



000; ^ (2Z + 1)! 

and then we substitute the Gosper approximation for all these factorials, that is, 



(44) 



/!~^(2/+i)^/'e-'. (45) 
The discrete multipole weight function (43) then reduces to a straightforward continuum version 



, 1 (2/i + l)(2;2 + l)(2/3 + l)(2;+|) {21 - 2h + l){2l - 2h + \){2l - 2h + \) 

w[li,l2,t3) (2/ _ 2h + i)(2Z - 2/2 + l)i2l - 2h + ^) V {21 + i) ^^^^ 

This is a remarkably accurate representation for the exact discrete wij^j^ with the difference between the 
weight functions being less than 0.01% (0.1%) for about 95% (99%) of the allowed triples 11^2^3 on the domain 
(42) with 2 < li, I2, 13 < 2000. The worst approximation by w{li,l2, 13) never differs by more than 2.5% and 
such points are exclusively located very near the boundaries, leaving an overall integrated error over the 
entire domain (42) of less than 0.01%. Nevertheless, care must be exercised using this approximation for 
edge- or corner-weighted models. With this caveat in mind, we can define the multipole sum equivalent to 
the wavenumber tetrapyd integration (41) as 

nf]= wimJim, = I f w{hMM)f{hMM)dVr, (47) 

with the inner product again defined by (/, g) = T[fg\. It will be clear from the context whether we are 
dealing with multipole or wavenumber integrations. 

The weight function w{li,l2,h) (or wi^i^i^) in (46) possesses an overall scaling which grows linearly with 
I, as illustrated in fig. 4. It can be convenient to eliminate this scahng, so that the weight function becomes 
very nearly constant. We can achieve this by dividing w{li,l2, h) by a separable function as 

n 1 1 \ w{li,l2,h) 
^^^^^'^^'^'^ -(2/1 + 1)1/3(2/, + l)i/3(2Z3 + l)V3- (^8) 

The result is shown in fig. 4 where it is evident that Ws ~ const, everywhere except on the boundaries. 
For uniform or centre-weighted bispectrum models, such as the equilateral model, the multipole domain 
with weight Ws{li,l2,l3) becomes essentially identical to that for the primordial wavenumbers (40) with 
w{ki, k2, = 1, so that it is a good approximation to proceed with the same polynomial expansions. 

Finally, we comment on the freedom to absorb an arbitrary separable function vi into the weight functions 
w{ki, k2, ks) or wi^i^i^, such as in the example (48) above. If we define a new weight w in the estimator as 

wiM, = Whi,iJ {vhVi^vi^f , (49) 



then we must similarly rescale the estimator functions as {bij^i^/A) = vi^vi^vi^{hi^i^i^ / A.) . This rescaling 
should be separable, otherwise it would compromise the separability of the methods we outline here, under- 
mining their efhciency. As we have seen it can prove convenient to make the weight functions scale-invariant 
for practical purposes, thus facilitating better convergence of mode expansions for typical bispectra. How- 
ever, in principle, we can also exploit this separability in order to remove pathologies from singular shapes, 
such as the local model, using a mode expansion to describe the more regular deviations away from these 
shapes. The important point is to consistently use both the new weight w and the estimator rescaling 
throughout the analysis pipeline, including the generation of appropriate orthonormal mode functions. 
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Figure 4: Scaling comparison of the multipole domain weight function w{li,l2,h) (or wi^^i^i^) given in (46) and the modified 
weight function Wa{li, h, h) given in (48), which is rescaled by a separable function. On the left, the equal-Z values are shown 
with the linear scaling of w (dashed) contrasting with the constant Wa (solid). On the right, a density plot of Ws is shown on 
the li + h + Is = 2L slice with L — 2000. Note the uniformity Ws ~ const., except very close to the edges where there is about 
a factor of 4 rise to the maximum value on the perimeter. 



B. Orthogonal polynomials on a tetrahedral domain 



We next construct some concrete realizations of mode functions which are orthogonal on the tetrahedral 
domain Vj- and which have the form required for a separable expansion (39). First, we will generate 
one-dimensional orthogonal polynomials qp{x) for unit weight w = 1, before discussing their promotion to 
three-dimensions and alternative weights. These tetrahedral polynomials are analogues of the more familiar 
Legendre polynomials Pnix) on the unit interval. Considering functions qp{x) depending only on the x- 
coordinate, we integrate over the y- and z-directions to yield the reduced weight function w(x) for x £ [0, 1] 
(we take K = 1): 



w{x) 



ix(4 - 3x) , with r[f] = C /( 

JO 



x) w{x) dx . 



(50) 



This simplifies our domain integration (41) for functions of only x, and the moments for each power of x 
become simply 



r[x-] 



2{n + Z){n + 2) 

From these we can create orthogonal polynomials using the generating function, 



(51) 



qn{x) 



1 

M 



1/2 7/24 1/5 
7/24 1/5 3/20 

Wn-l Wn Wn+l 



1 



X 



X 



Wn 
Wn+l 

W2n-1 



(52) 



where we choose the normalisation factor M such that T[q„ 
are orthonormal 



1 for all n G N, that is, so that the qn{x) 

'^[q-nqp] = / qn{x) qp{x) dVr = Kp ■ (53) 

JVt 
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Figure 5: The orthonormal one-dimensional tetrahedral qn(x) plotted on the unit interval for n = 0-5. The behaviour is smooth 
and bounded across the domain even for high n, except where the weight function w{x) vanishes at a:: = 0. Also plotted for 
comparison are the rescaled Legendre polynomials Pn{2x — 1) (dashed lines). Despite qn and P„ sharing qualitative features 
such as n nodal points, their properties and orthogonality on Vt are very different. 



The first few orthonormal polynomials on the tetrahedral domain (40) are explicitly 

qo{x) = V2, 

qi{x) = 5.787(-^ + x), 

q2{x) = 23.32 fx + , (54) 

q^{x) = 93.83 (-0.09337 + 0.7642 X- 1.631 x^ + x^) , 

qA{x) = 376.9 (0.03192 -0.4126X + 1.531x2 -2. 139 x^ + x"^) , 

q^ix) = 1512 (-0.01033 + 0.1929 X- 1.084x2 + 2.549 x^- 2.644 x^ + x^) , ... 

These can be obtained easily from the generating determinant (52) in Mathematica or similar applications. 

We note that the g„'s are only orthogonal in one dimension (e.g. we have T[g„(x) qp{y)] 7^ 5np in general). 
However, as product functions of x, y and z they form an independent and well-behaved basis which we 
will use to construct orthonormal three-dimensional eigenfunctions. In practice, these qnS will remain the 
primary calculation tools throughout, notably when performing separable integrations. Where they differ 
from the separable functions used to represent bispectra in the literature, they generally have a number of 
distinct advantages, as we shall detail at the end of this section. Finally, we point out that for a regular 
tetrahedron (in contrast to the tetrapyd domain (40)), the volume weight function is tD(x) = 2x(l — x) 
and so the behaviour is different at x = 1 where the weight vanishes, unlike (50). The first orthonormal 
polynomials in this case are go (x) = ^/3, qi{x) = 0.387(2y - 1), 52 = 32.4(y - 0.724)(y - 0.276) , ... . 

ow let us turn to the polynomials q{x) which are orthonormal on the multipole domain (42), using the 
weight functions w given in (46) and Wg given in (48). For definiteness we take L = /max = 2000, so that 
X = li/L, y = I2/L and z = l^/L. The generating function (52) can be obtained as above but now using 
the moments Wn = T[x"'] = J w{x, y, z) x^'dVq- (or by undertaking Gram-Schmidt orthogonalisation from 
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Figure 6: Orthonormal polynomials qn{x) for the multipole domain (42) with weight functions w given in (46) [solid] and Ws 
given in (48) [dashed], as well as the previous qn(x) for unit weight [dashed] (shown already in fig. 5). Despite the different 
scaling of w, these tetrahedral polynomials are very similar, particularly the latter two with flattened weight functions. 



qq = const.). The resulting first few polynomials for the multipole domain are then 





= 0.07378, 




qi{x) 


= 0.3017 (-0.6110 + a;), 




q2{x) 


= 1.223 (0.2665 - 1.145 x + x^) , 






= 4.933 (-0.1000 + 0.7951 X- 1.659x2 


+ x') , 


Q4{x) 


= 19.85 (0.0345 - 0.4342 x + 1.578x2- 


2.169x3 + x'^) , 


Q5{x) 


= 79.55 (-0.0106 + 0.1975 X- 1.103x2 


+ 2.576x3 -2.657 x^' + x 



A cursory comparison with qn given above for the flat wavenumber domain will show that these polynomials 
are very similar for low n, despite the linear scaling behaviour of w. However, if we remove this scaling as 
in the flatter weight Ws in (48s), the polynomials become near identical as illustrated in fig. 6. It is clear 
that each of these polynomial sets would suffice as independent basis functions on the multipole domain. 
However, using the correctly weighted versions leads to improvements in the immediate orthogonality of 
the three-dimensional polynomials we shall construct in the following discussion. 



C. Bispectrum symmetries and three-dimensional basis functions 

We can represent arbitrary bispectra on the tetrahedral domain (40) using a suitable set of independent 
basis functions formed from products qp{x) qr{y) Qs{z) of the orthogonal polynomials (54) (or with different 
weight functions, such as (55). (Here, we again take x = ki/kj^ax, V = ^2/fcmax5 z = ks/k^riax or x = Ai/Zmax, 
etc.) Both primordial bispectra B(ki, k2, k^) and CIVIB bispectra bi^i^i^ on (40) possess six symmetries made 
from combinations of discrete tt/3 rotations around the line x = y = z and/or reflections which interchange 
the axes. We can impose these six symmetries on our products by summing the relevant permutations and 
defining the 3D basis function 

Qn{x,y,z) = [qp{x)qr{y)qs{z) + qr{x)qs{y)qp{z) + qs{x)qp{y)qr{z) 
+ qp{x)qs{y)qr{z) + qs{x)qr{y)qp{z) + qr{x)qp{y)qs{z)] 
= q{p qr qs} with n ^ {prs} , (56) 



16 



where we use the notation {prs} to denote the six permutations of pr.s. Here, for convenience, we have 
specified a one-to-one mapping n ^ {prs} ordering the permuted indices into a hst labelled by n (see below) . 
Alternatively, we could directly represent bispectra in a power series using sums of monomial symmetric 
polynomials which like (56) are also separable; that is, we could identify our set of basis functions with the 
following 

1, x + y + z, xy + yz + zx, x^ + y^ + z^, xyz , x^ + y^ + z^, etc. (57) 

The Qn{x, y, z) we defined in (56) are themselves ultimately constructed from these through the qp products. 

However, the Q„ have two distinct advantages which are, first, they already have partial orthogonality 
built in which improves their convenience and convergence and, secondly, unlike the elements of (57), 
the qp polynomials remain bounded and well-behaved when convolved with transfer functions, as we shall 
emphasise in the map-making discussion. 

Since we will be dealing with relatively small numbers of basis functions, it is convenient to order the 
symmetric products Qn = q{p qr qs} linearly with a single index n; here we offer two comparable alternatives 
for achieving this. The first is by 'slicing' such that triples are ordered by the sum p + r + s and the second 
is by 'distance' from the origin, that is, p'^ + r'^ + s'^. 

Slicing the prs naturally groups the Qn by the overall order of the polynomials from which they are made. 
The subscript n, with a specific choice of sub-ordering, relates to the prs via 



^ 000 


4 ^ 


111 


8 


022 


12 ^ 


113 


1 001 


5 ^ 


012 


9 


^ 013 


13 


023 


2 ^ Oil 


6 ^ 


003 


10 


^004 


14^ 


014 


3^002 


7^ 


112 


11 


^ 122 


15^ 


005 



where we have underlined the transitions between polynomial order. The number djv of independent 
symmetric polynomial products QnQpQr which can be formed at each polynomial order N isa, combinatorial 
problem but the sequence begins as follows and we give a recurrence relation for any further elements: 

{d7v} = {l,l,2,3,4,5,7,8,10,12,...}, djv = 1 + djv-2 + - ^JV-s ■ (59) 

For consistency when using slicing we will usually decompose functions with polynomials up to a specific 
order A'^. 

The distance ordering of the Qn is more straightforward with 

0^000 2^011 4^002 6 ^112 8^122 

1^001 3^111 5^012 7^ 022 9 ^ 003 • • • . (60) 

This approach is the analogue of state counting over spherical shells in the continuum limit and the basis 
functions can be grouped accordingly. Distance ordering has some advantage by reshuffling to higher n the 
pure states 00^» which turn out to be most affected by masking. 

While the Q„'s by construction are an independent set of three-dimensional basis functions on the domain 
(40), they are not in general orthogonal. In fig. 7, we illustrate the inner product matrix = (Q„, Qp), 
showing partial orthogonality (nearly diagonal 7„p) because of their origin as products of orthogonal q^s. 
However, this is not sufficient because we need the convenience of a fully orthonormal basis to efficiently 
decompose arbitrary bispectra. For this reason, we undertake an iterative Gram-Schmidt orthogonalisation 
process to construct an orthonormal set TZn from the Qn, that is, satisfying 

{Un,Up) = 6np. (61) 

Formally, we have a Gram matrix T = {{Qn, Qp)) made from the independent functions Q„, and therefore 
positive definite, which needs to be factorised as F = A'''A where A = {{Qn, T^p)) is triangular (i.e. an LU 
or Cholesky decomposition). As we require explicit relationships between Qn and TZn, we run through the 
main steps in the Gram-Schmidt process. 
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Figure 7: Partial orthogonality of the symmetric product polynomials Q„ illustrated through the inner product matrix {Qn, Qp) 
for < n,p < 10 (left panel). Lower triangular matrix A„p in (62) illustrating the decomposition of the orthonormal TZn into 
the Qp arising through the Gram-Schmidt process (right panel); this is the inverse of {Qn, TZp). To improve comparison, the 
Qn's have been unit normalised. 



Let us assume that we have achieved this orthonormaHsation up to n, that is, such that {TZn, Tim) = Snm-, 
y m <n. This means we can represent any basis function Qp in terms of the IZm and vice versa by inversion, 
so we can write 



T^m = ^ AmpQp for m,p<n, 

p=0 



(62) 



where Xmp is a lower triangular matrix with (A 



Inp 



{Qn-, Tip)- We wish by induction to construct the 



next orthonormal polynomial IZn+i and infer from this the sum over basis functions up to Qn+i- We achieve 
this by taking the next independent basis function, Qn+i, as a first approximation to an unnormalised "R-'n+i 
and then we project out all components dependent on the TZm {m < n), 



n+l 



IL-T-L It' n 

T^ri+l = ^A^+ipQp = Qn+1 - X] Tim I Qn+lTimW dVr 
p=0 m=0 -^^^ 

n m m 

= Qn+1 — Xmr'^ms^n+l s Qr 



(63) 



m,=0 r=0 s=0 



where in the second line we have substituted (62) and the 7n+is are determined from the relative orthogo- 
nality of the Qn's, 



In+ls = (Qn+1, Qs) = j Qn+1 QsWdVr 

By equating coefficients in the expression (63) we can determine that 



(64) 



A'ra+l p — ^n+lp — ^ ^ Kp Ks 7n+l ; 
r=p s=0 



(65) 




Figure 8: Three-dimensional orthonormal polynomials TZn on the tetrahedral domain (40). Taken from top left (and moving 
across and then down) these are TZo, TZi, TZ2, TZs, TZa, and TZai (bottom right). 
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Unit normalising appropriately, we obtain the coefficients A^+ip which define the new orthonormal TZn+i 
we are seeking, that is, we have 




In fig. 7, we see the orthogonalisation process at work for the first 10 modes by plotting the matrix 
coefficients for {Qn, Qp) and (Q„, TZp). At each order n, the independent component in TZn is provided 
by Qn, as indicated by the dominant diagonal term. This is a good approximation at low order, but the 
mixing increases with n. We also illustrate several of the orthogonal polynomials TZn on the tetrapyd 
domain fig. 8 for the slicing ordering (58). These are primarily the lowest modes and demonstrate the build 
up of the number of nodal points and lines as the order increases. As an aside, we note Gram-Schmidt 
orthogonalisation in the form given above is inherently unstable numerically, though this can be easily 
corrected by using the modified Gram-Schmidt process. However, we do not iterate to sufficiently high n 
to notice any significant degradation in accuracy, as verified by determining orthogonality. 



D. Mode decomposition of the bispectrum 

We have constructed examples of an orthonormal basis {TZn} out of monomial symmetric polynomials 
(57) which span the set of symmetric functions on the tetrahedral domain (40). The TZn polynomials will 
possess the properties of more familiar orthonormal eigenmodes in other contexts, notably completeness and 
the convergence of mode expansions for well-behaved functions. We proceed by considering an arbitrary 
primoridal bispectrum (12) described by the shape function S{ki, k2, ks) and decomposing it as follows 

oo 

S{ki,k2,ks) = '^anTZn{x,y,z) , (67) 

n=0 

where the expansion coefficients are given by 

< = (7^„, 5) = / TZnSw dVr , (68) 

JVt 

and K = fcmax and ki = Kx etc on the domain Vr defined in (40). For practical purposes, we shall always 
work with partial sums up to a given N = rimax with 



N 

Sn = y^a'^TZn{x,y,z) , S= lim Sn ■ (69) 

AT— »oo 



n=0 



We shall assume that the expansion (69) is the best fit mode expansion of degree N (for this particular 
mode ordering). Given the complete orthonormal basis TZn, Parseval's theorem for the integrated product 
of two functions implies 



N 

{S,S')= / SS'wdVT= hm (70) 

/li N—*oo ' 

n=0 



which, for the square of a function S, yields the sum of the squares of the expansion coefficients, Tf/S^] = 

Z^n n 

In order to accomplish our original goal of a general separable expansion (39), we must now transform 
backwards from the orthonormal TZn sum (69) into an expansion over the separable product functions 
Qn = qipQrQs} through 

N 

Sn = J2^nQn{x,y,z) , (71) 

n=0 
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Figure 9: Correlation of the reconstructed bispectra to the original for partial sums of the decomposition up to a given mode n. 
The plot includes the primordial bispectra for the equilateral and DBI models, the CMB bispectrum for the equilateral and DBI 
models and the CMB bispectrum produced at late times by cosmic strings. In all cases, we find that with 15 three-dimensional 
modes we have a correlation greater than 98%, thus demonstrating very rapid convergence. For the CMB bispectra, convergence 
is limited by matching the acoustic peaks introduced by the transfer functions, whereas the primordial models converge at 98% 
accuracy with only 6 modes. 

where the a® can be obtained from the as 

N 

a^ = ^{\')npa;, (72) 
p=0 

with the transformation matrix Xnp defined in (62) (this is triangular and not orthogonal in general). Note 
the complication that also contains contributions from TZp components with n < p < N, since {X^)np is 
upper triangular. The inverse transformation 

N 

«n =E(^"')np«p. (73) 
P 

has coefficients given by (A~^)„p = (Q„, TZn)- We have already noted that the degree of non-orthogonality 
of the Qn basis is described by jnp = {Qm Qp) in (64) which is in turn related to Xnp through 

N 

(7"')np = J^(A^)nrA.p . (74) 

r 

When substituted into Parseval's theorem (70) in the Qn basis, we see that the coefficients of different 
degrees become mixed as 

N N N 

{Sn, S'n) = XI ""^ = X] 2Z "nlnpOp (75) 
n n p 

The separable Qn expansion (71) is important for most practical calculational purposes but its coefficients 
are constructed at the outset using the orthonormal TZn- For interpreting results from the estimator it 
is helpful to transform back to the TZn basis in order to understand the normalised spectrum using 
Parseval's theorem (70). We finally note that all the transformation matrices, Xnp and jnp in (64), need 
only be calculated once, at the same time as the TZn polynomials are generated, and then stored for later 
reference. 
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Figure 10: Decomposition into orthonormal polynomials TZn for both the primordial shape function (67) (above) and the CMB 
bispectrum estimator (89) (below) for the equilateral (red) and DBI (blue) models. In both cases, these results are for 'slicing' 
polynomial ordering given in (58). The peak in the CMB bispectrum estimator modes (here, for /max = 500 at n ~ 5) arises 
because of the power shifted into the coherent acoustic peaks observed in fig. 3; this is a distinguishing feature of the CMB 
bispectrum for most primordial models [14]. 



In fig. 9, we demonstrate polynomial convergence for the DBI model and its separable equilateral approx- 
imation by showing the cross-correlation between the shape function and the partial sum (69). We also 
provide the actual expansion coefficients for the primordial shape functions in fig. 10 (along with the 
for the CMB bispectra). Using only 6 three-dimensional IZn polynomials we achieve a better than 98% 
cross-correlation with the original analytic expressions in both cases (i.e. using symmetric products of at 
most quadratic qp polynomials from (54)). Here, we undertake the full Fisher matrix analysis between the 
theoretical CMB bispectrum and its approximation using the methods described in ref. [11]. More generally, 
we note that for all well-behaved bispectra the polynomial expansion has proved to be rapidly convergent. 
The decomposition of the primordial bispectra is also numericaly efficient using the orthogonal IZn modes 
with each coefficient taking an average of 7 seconds to calculate. 

We can equally well expand the CMB bispectrum h-^i^i^ at late times, using the same polynomials h-^i^i^ = 
a^7^„(x, y, z). However, our aim is to represent the bispectrum estimator £ given in (102), rather than 
6(^2/3 itself. We, therefore, consider expanding a separable product which approximates £ with the same 

weight and scaling (schematically, \/l bi^i^i,JCf^'^). We discuss this in the next section, but in the lower 
half of fig. 10 we show the corresponding late-time expansion coefficients for the equilateral and DBI 
CMB bispectra. Once again, convergence is rapid, see figure 9, with a 95% correlation achieved with only 
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12 TZn polynomials at Inrnx = 500, improving to 98% with 15 polynomials for both CMB bispcctra. This is 
despite the fact that the expansion must incorporate additional features induced by the transfer functions. 
We emphasise the power shift from the low modes in the primordial bispectrum to a peak at higher modes 
n ~ 5 in the CMB bispectrum (for this slicing and /max)- This is a common characteristic of the polynomial 
expansion for almost all bispectra of primordial origin and is a manifestation of the pattern of coherent 
acoustic peaks observed in ref. [14]. The decomposition of the CMB bispectra is also numericaly efficient 
using the orthogonal TZn modes with each coefficient taking an average of 8 seconds to calculate. 

E. Utility of the tetrahedral polynomials Q„ and other alternatives 

The three-dimensional polynomials Qn we have presented are just one possible set of basis functions which 
can be used as bispectrum eigenmodes for the methodology we present in the next section. They are built 
from products of the one-dimensional g^'s which are orthonormal on the tetrahedral region (40) with given 
weight functions. These are the analogues of Legendre polynomials Unfortunately, unlike the P„'s on a 
cube, they do not retain full orthogonality as separable products on the tetrahedral domain, though there is a 
substantial remnant. There are significant advantages to using the Qp's, rather than the monomial symmetric 
polynomials in (57), in the same way that Lcgcndrc polynomials arc more efficient than simple power scries 
representations. As we shall discuss subsequently, there are further important benefits which arise when 
the Qn's are decomposed into separable integrals over the Qp's. Given the bounded and well-behaved nature 
of the Qp's on their domain, these integrals reflect these properties, eliminating diverging artifacts which 
are known for other separable approximations to bispectra in the literature (including difficulties for simple 
powers x"-). 

There are other alternatives to expansions using the tetrahedral polynomials Q„ and TZn which we have 
considered. It is possible, for example, to expand an arbitrary bispectrum using separable products of more 
familiar orthonormal functions such as Legendre the P„ and Chebyshev Tn polynomials, as discussed in 
ref. [14]. This entails using shifted polynomials on the full cubic domain li,l2,h < ^max- The shortcoming 
of this approach is that the bispectrum is only defined on the tetrahedral region (40) , so it has to be zero 
elsewhere or arbitrarily extended in some manner to fill the cube. This leads to generic overshooting of 
the expansion near the boundaries (the analogue of the Gibbs phenomena for Fourier series). Extensive 
experiments yielded very poor convergence with Legendre and Chebyshev polynomial expansions, as well 
as Fourier series, especially relative to that achieved with the tetrahedral Q„ and TZn polynomials. A 
further simple alternative is to transform the tetrahedral region into a cube (see ref. [11]). This allows 
the bispectrum to be defined everywhere on the standard domain using the more familiar eigenmodes and 
thus yielding more rapid convergence. However, this compromises separability which is essential for the 
estimators we discuss below. 

Were the rate of convergence to become a primary issue when representing the bispectrum, then there 
are further alternatives to polynomials. There is a significant literature on eigenmodes on the regular 
tetrahedron or simplex because of its importance in crystallography and other contexts. For example, 
it is possible to define generalised sine and cosine functions on the simplex, as well as Koornwilder and 
generalised Chebyshev polynomials of the first and second kind (see, for example, ref. [28]). Such generalised 
eigenfunctions could, in principle, improve convergence, however, two significant developments are required. 
First, it is more natural to define the observational data on the tetrahedral domain with li,l2, h < ^max (the 
tetrapyd), rather than the simplex I1+I2+I3 < 2/inax; so generalised eigenfunctions must be derived explicitly 
for this domain (42). Secondly, these should be able to conveniently represent functions in separable form. 
The present tetrahedral polynomials TZn and Qn do converge satisfactorily for all the primordial models 
studied to date, but more efficient mode expansions will continue to be investigated [29]. 
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IV. MEASURES OF Fnl 
A. Primordial Fnl estimator 

We have obtained two related mode expansions for a general primordial shape function (12), one for an 
orthonormal basis TZn (69) and the other for separable basis functions Qn (71). Substitution of the separable 
form into the expression for the reduced bispectrum (14) offers an efficient route to its direct calculation 
through 



1I2I3 



(|)^A|/jvL / x''dxdkidk2dk3 6J2< Qn(A;i,A;2,A;3)Ai,(/ci)A,,(A;2)A,3(A;3)ii,(A;ix)ii,(A;2x)ii3(fc3x) 
x'^dxlf^ J dki qp{ki) Ai^{ki) ji^{kix)^ ijr J ^^^^ q'r-(fc2) Ai2(A;2) ji2(A;2x) 



J dksqsiks) Ai^{k3) ji^{k3x)j + 5 permutations 
= A|/jvL^aS J x^dxq\^/,^ql^ = A|/^i^a« j x^dx 2^^'% 



(76) 



where here we implicitly assume the mapping n <-> prs between indices for the Qn and the product basis 
functions from which they are formed, that is, Q„ = q{pqrqs} (e-g- see the ordering in (58)). For brevity 
we have also denoted as the convolution of the basis function qp{k) with the transfer functions 

(77) 



qj,{x) =-j dkqpik) Ai{k)ji{kx) , with Qll^^^'ix) = q[l{x) ql^{x) q^^fix) . 

(These gr^ are the primordial counterparts of the q defined in multipole space (55).) Here, in (76), the previ- 
ously intractable three-dimensional wavenumber integral separates into the product of three one-dimensional 
integrals which are relatively easy to to evaluate. This has been achieved because the triangle condition has 
been enforced through the product of Bessel functions, giving a manifestly separable form and allowing us 
to interchange the orders of integration with x; it is the basis for the analytic local (19) and constant (16) 
solutions on large angles, as well as all the analysis of separable shape functions to date (see, for example, 
ref. [9]). With this mode expansion, all non-separable theoretical CMB bispectra bij^^i^ become calculable 
provided there is a convergent expansion for the shape function. Accurate hierarchical schemes already 
exist against which to benchmark this method [14] but, in principle, it is more efficient. 

Now consider the implications of this mode expansion for /nl by substituting the decomposed bi^j^ (71) 
into the estimator expression (31) to obtain 

S = Alfr,Lj2<Il I x'dxql;^{x)qHx)qil{x) J d'n Y^,^, (n) 1^,^, (n) Y,,^, (n) "'^'g »^-^^^3n.3^ ^8) 



li,mi 



li,mi 



Q2 



l3,m3 



CI3 



(79) 



The break up of the wavenumber integration now extends also to the separation of the sum over the 
multipoles h,l2,h- The summation between the aim's and each qp integral creates a filtered map of the 
original data, which we can define in the above as 



Mp(n,x) = ^g. 



m 



TT 



qp{k)Ai{k)ji{kx)dk 



aimyimjii-) 

Ci 



(80) 



Im Im 

From these we can efficiently calculate product maps which essentially extract the Qn basis function con- 
tribution from the observational data, 

M^{n, x) = Mp(n, x)Mr{h, x)Ms{n, x) (81) 
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where again we exploit the correspondence n prs. Note, that in this case, there is no need to symmetrise 
the product map because it is imphcit in the estimator expression. Integrating over directions and shells 
we can now obtain for the observational maps, the analogue of the primordial mode expansion coefficients 



a® 



(fh / x'^dxM^{h,x). (82) 



In common with the analysis of simple separable shapes, the shell integral over x is where the most significant 
computational effort is required. 
Substituting into (78), the bispectrum estimator then collapses into a compact diagonal form 

S = ^J2<^n- (83) 

n 

The estimator has been reduced entirely to tractable integrals and sums which can be performed rapidly 
even at Zmax = 2000. We will demonstrate how efficiently it can recover /nl from simulated maps in 
subsequent sections. 

The form of the estimator in (83) suggests that further information can be extracted about the observed 
bispectrum beyond the /nl for one specific theoretical model. This is because, through the coefficients 
we have obtained some sort of mode decomposition of the bispectrum of the observational map. However, 
the non-orthogonal and primordial nature of the Q„ basis functions means that these require some 
effort in their interpretation. Consider the expectation value of /?® obtained from an ensemble of maps 
generated for a particular theoretical model with shape function S = J2n'^n Qn- Noting that the relation 
miO-i2rn2^lsm3) — ^mim^mj^hhh-i average over the product maps (81) becomes 

m = Jd'hJ x'dx{MS{n,x)) = E (/ ^'dxq'l/^^q';^^ (^gl^J^ls^^y ^^^^ (84) 

111213 p 

where we have substituted the expression (76) for the reduced bispectrum and the weight wij^i^ is described 
in (43). Here, the matrix Tnp represents a late time inner product Tnp = {{Qn, Qp)) analogous to 7„p = 
(Qn, Qp) in (64) (but with a different weight so that ((7?.„, T^pJi 7^ ^np)- Determining the transformation 
matrix F.^^ relating the af^ and appears to be a complicated task but, in fact, it reduces to separable 
sums and integrals over the one-dimensional products qpQr convolved with Bessel and transfer functions, 
together with the final sum over the multipole domain (42) . The latter is straightforward, especially in the 
continuum limit (46). It need only be evaluated once, given a robust prior estimate for the power spectrum 
Q's. 

This discussion demonstrates that we can recover spectral information about the primordial shape function 
from the observational data through the relation 



= = ^(^r2-^ ) (/?«), (87) 
p 



np 



which extends to the orthonormal coefficients a'^ using (73). If the decomposition coefficients /3® are found 
with adequate significance, we can reconstruct the shape function from a single realization through the 
expansion 

Siki,k2,h) = J2{^^~')^ PpQn- (88) 
n,p 
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Wc will discuss this prospect in more detail in the next section about the late-time CMB estimator where 
the relation between the a® and /3® is more transparent. 



B. CMB Fnl estimator 

Wc turn now to the implementation details of the late-time CMB estimator originally proposed in [14]. 
Here, wc presume that the CMB bispectrum bi^i^i^ for our non-separable primordial model is precomputed 
using eqn (76) or a robust hierarchical scheme [11, 14]. In addition, this approach can accommodate any 
late-time source of non-Gaussianity in the CMB, including secondary anisotropies, gravitational Icnsing, 
active models such as cosmic strings, and even systematic experimental effects. For the late-time analysis 
we wish to expand the estimator functions using the orthonormal TZnihjhjh) and separable Qn(^i)^2)^3) 
mode functions created out of products of the qp{l) polynomials, for which we gave a concrete example (55). 
(Note that we denote the multipolc modes with a bar, distinguishing them from the primordial qp, Qn, TZn 
which are functions of wavenumber k). Convergence of mode expansions on the multipole domain (42) has 
been found to be poor for quantities as scale-dependent as bi-^i^i^, so we choose to decompose the estimator 
functions directly as 



biMs = 2^ot^Qn, (89) 

n 



where the separable vi incorporates the freedom to make the weight function wi^i^i^ given in (43) even more 
scale invariant (typically we shall use vi = {21 + 1)1/6 as defined in (48)). The expression (89) means that 
we are effectively expanding in mode functions modulated by the Q's, that is, Qn — > y/CiQn/vi- These 
more closely mimic the acoustic peaks observed in the bij^i^ as illustrated in fig. 11. We shall see that the 
estimator expansion with Ci in (89) is appropriate for primordial models, but different flatter choices will 
be more suitable for late-time anisotropy, such as that from cosmic strings. 

We determine the implications for /nl of our mode expansion (89) by substituting into the estimator 
(31), 



li,mi n'r-*prs 



(90) 



n^prs 'i V 'i y \l2,m2 V '2 y \h,m3 V J 



where again we assume the correspondence between the label n and an ordered list of permuted triples 
{prs}, through Q„ = q^pQrqs}- As previously for the primordial estimator (78), we note that the sum 
between the qp{l) and the aim creates filtered versions of the original CMB map defined by 

Mp(n) = ^q^(l)^^Yim{h) , (92) 

which are multiplied together in (90) to form the product map 

Mn{n) = Mp{h)Mr{n)Ms{n) . (93) 
Integrating over directions, we can obtain the map mode expansion coefficient 

Pn = J cPnMnin) . (94) 

Thus the estimator reduces again to diagonal form 

^=nT. <~Pn ■ (95) 

n=0 




Figure 11: Polynomials TZn on the tetrahedral domain (42) used for representing modes in the CMB estimator multiplied by 
the weight function given in (89). These are ordered just as in fig. 8 with IZq (top left), IZi, 1^2, TZ-i, IZ4,, and 7^.41 (bottom 
right). The last higher mode bears a superficial resemblance to the equilateral bispectrum in fig. 3. 
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Like (83), it consists entirely of separable sums and tractable integrals which can be performed rapidly. 

As before, the separation of the estimator into two and /3® halves indicates that this approach could 
offer more direct information about the bispectrum than just /nl for one model. Consider the expectation 
value of (3^ from an ensemble of maps with a given bij^h (and Ci) expanded as (89). Following the steps 
used to derive (84), we find a considerably simpler expression after substituting (89): 

= E Y^^mM)^^^^^^ or) 

C;:l^vi,vi,vi,^Ci,Ci,Ci, y ^ Vhvi^vi, 

P lihh P 

where the modified weight function wi^i^i^ is given in (49) and T^p = {Qn, Qp) as discussed previously. 
Hence, the estimator, when applied to a map containing the bispectrum defined by a®, should have the 
expectation value 

(^) = ^EE"nfnp«=. (99) 
n p 

Now rotating to our orthonomal basis we note that from the relation (??) we can deduce the simple 
and elegant form 



AT 

n 

That is, we expect the best fit P^^s for a particular realization to be the a^'s themselves. The simplicity 
of this result is not unexpected, since it would be obtained by correlating a bispectrum decomposed into 
the TZn with itself. The advance here is that extracting the spectrum from the observed map would be 
intractable for large /max; were it not for the transformation made to a non-orthogonal separable frame. 
Assuming the coefficients are measured with some significance from a particular experiment, we can go 
further and reconstruct the map bispectrum using (89) 

We reiterate that the viability of this fast and general reconstruction scheme [14] depends on two key 
factors, first, the smoothness of the reduced bispectrum hi^i^i^, requiring few modes to characterise it, and, 
secondly, on the completeness of the orthonornial basis from which the separable expansion was obtained. 
We note that this methodology can be applied using any complete mode expansions, beyond the polynomial 
examples given here, as well as with over-complete decompositions, such as wavelets, or with binning. In 
the next section, we will demonstrate the efficacy of this method with simulated maps (for a sufficiently 
large /nl), recovering the expected spectrum and the main distinguishing features of the bispectrum 



C. Observable Fnl normalization 



In previous work [11], we pointed out the shortcomings of normalising the quantity /nl using the conven- 
tions employed to date in the literature (see also [25]). At present, the central point in the primordial shape 
function defined in (12) is normalised to unity assuming scale invariance, that is, S{k,k,k) = 1 with no 
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^-dependence. This produces inconsistent results between models peaking or dipping at this central point 
(actually along this line); contrast the factor of 7 between the quoted variances of the equilateral and local 
models for exactly this reason. Furthermore, the definition is not well-defined for models which are not 
scale-invariant, such as feature models, and it is simply not applicable nonGaussian signals created at late 
times, such as those induced by cosmic strings or secondary anisotropics. 

We, therefore, propose a universally defined bispectrum non-Gaussianity parameter Fnl which (i) is a 
measure of the total observational signal expected for the bispectrum of the model in question and (ii) 
is normalised for direct comparison with the canonical local model (in particular, with F^f^ = f^f^ for a 
given imax)- We presume that we have an unnormalised CMB bispectrum bi^^i^i^ accurately calculated for 
a specific theoretical model over the whole observationally relevant domain I < /max- This can be achieved 
for any model using the separable mode expansion (76) or hierarchical methods [11]. We then define -Fnl 
from an adapted version of the estimator (31) with 

where N is the appropriate normalisation factor for the given model, 
and iVioc is the normalisation for the local model with /nl = 1, 

^loc (/nl=1)^ 

Nil = E r ■ (1°^) 

This Fnl estimator will certainly recover the usual /nl for the local model, but it is also clear that it will 
also equitably compare the total integrated observational bispectrum with that obtained from the /nl = 1 
local model. Of course, these definitions presume a sum to a given / = /max (which should be quoted) 
but results for primordial models should not depend strongly on this cut-off, unless scale-invariance is 
broken. In any case, diffusion from the transfer functions means that the primordial signal is dying out 
beyond / > 2000, so we propose a canonical cut-off at /max = 2000 (which is also relevant in the medium 
term for the Planck experiment). Late-time anisotropics, such as cosmic strings, do not generically fall-off 
exponentially for / > 2000, but meaningful comparisons to the local /nl=1 niodel can be made with the 
same definition (102) on this domain, and alternative measures can be proposed elsewhere. In principle, the 
normalised estimator (102) can also be adapted as a gross measure of the total bispectral signal over the 
given domain, irrespective of the possible underlying physical model. For example, using the reconstruction 
from Parseval's theorem (70), the estimator provides a measure of which should then be normalised 
relative to the total expectation for the local model with N = N\oc in (102). 

If the CMB bispectrum h^i^i^ is not known precisely for the primordial model under study, then the 
normalisation factor N in (103) can still be estimated using the shape function S{ki,k2jkz)- Primordial 
and CMB correlators are closely related, so one can obtain a fairly accurate approximation to the relative 
normalisations above (103-104) from [11] 



= / S'^{ki,k2,k^)w{ki,k2,k^)dVk, (105) 



where the appropriate weight function was found to be ■w{ki,k2,k^) ~ l/(^i + k2 + k^) and the domain 
Vk is given by ki,k2,k3 < A;max(/max) (refer to the discussion before (28) in section II). Here, we note 
that A^/iV;oc/NL=i ^ N /Niocf^^=i- Using this primordial shape function normalisation N in ref. [11] led 
to a comparable definition of /nl ~ Fnl, which can be useful for making fairly accurate projections 
of nonGaussianity or for renormalising /nl constraints for different models into more compatible Fnl 
constraints. 
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Below wc renormalisc published and model-dcpcndcnt constraints on /nl into the integrated bispectral 
measure Fnl, using the expression (102) together with accurate calculations of Bij^i^ for each model: 



-AKf^l'^^KSO [7] -4<Fi[fL < 80 (106) 

-125 < < 435 [10] ^ -24 < < 83 (107) 

-375 < f^l"'' < 37 [30] -93 < F^l'"' < 9 (108) 

-369 < f^'l^" < 71 [10] ^ -114 < F^'l^" < 22 (109) 



Note the much more consistent variance found for the different models with -Fnl, thus aiding direct com- 
parison, as well as the exact correspondence for the local model to which it is normalised. 

V. CMB MAP SIMULATIONS FOR GENERAL BISPECTRA 
A. Map-making with separable shape functions and its limitations 

In the limit of weak non-Gaussianity, an algorithm to produce non-Gaussian CMB simulations with a 
given power spectrum and bispectrum for separable primordial shapes was described in ref. [31]. Here, we 
present it in a more transparent notation, generalising the method to non-separable shapes using the mode 
decompositions of the previous sections. A byproduct is that the generalised approach is more robust and 
reliable, because the polynomial mode functions are better behaved than the separable approximations which 
have been previously employed. In this algorithm the non-Gaussian components of the CMB multipoles are 
obtained using the following formula: 

hrrii ^ / 2 3 

where is the Gaussian part of the CMB multipoles, generated using the angular power spectrum Q, 
while -8^2/3 is the given bispectrum of the theoretical model for which simulations are required. Although 
equation (110) is completely general, as before, its numerical evaluation is only computationally affordable 
for bispectra that can be written in separable form. We have emphasised already that separability results in 
a reduction of the computational cost of the estimator (31) from 0{l^^^) to 0{lf^g^^) operations; the same 
argument applies here allowing a rewriting of (110) into an equivalent form in pixel space (see below). 

The limitation dictated by separability is clearly overcome by using our cigcnfunction representation 
for the bispectrum (71). The basic idea is to start by expanding an arbitrary bispectrum shape S using 
the separable polynomial decomposition Sn until a good level of convergence is achieved and then to 
substitute the mode decomposition into (110). The accuracy of convergence is parametrized in terms of 
the correlation C{S, Sjsf) between the original non-separable shape and the eigenmode expansion, as defined 
previously (28). Note that this convergence can also be checked more accurately using the full Fisher 
matrix correlation on the CMB bispectra C{bij2i^,bi^i^iJ, calculated using the separable approach (76) or 
else accurate hierarchical approaches [11]. In previous sections (see fig. 9), we have noted how rapid this 
convergence is for well-behaved non-separable shapes, such as DBI inflation (or at late times with cosmic 
strings) . 

In addition to the bispectrum separability requirement, there is an important further caveat which can 
prevent the straightforward implementation of the algorithm (110). By construction, terms C(/nl) 
higher are not explicitly controlled. Following the discussion in [32] we can write the connected N-point 
functions as: 

(aU,«U) = K + fLc[^''] (111) 
{al^y^^^^a'^"^^) = [/NLi?/,M3+0(/NL)] (112) 

^^iimi^Z.m2^i3m3 _ _^/„m„^ = O (/^J . (113) 
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Thus the condition that the map has the power spectrum C; specified in the input will only be satisified 
if the power spectrum of the non-Gaussian component in (111) remains small. Since this method does not 
control 0(/nl) terms, one has to ascertain that spuriously large C/^*^ contributions do not affect the overall 
power spectrum significantly. It turns out that this effect plagues current map simulations if the standard 
separable expressions for the local and equilateral bispectra are directly substituted into (110), as we now 
demonstrate. 

In section II, we showed how the reduced bispectrum could be written explicitly in separable form for the 

local model in (20) and for the equilateral model b'^'j^f^ in (23). These were expressed in terms of one- 

dimensional convolution integrals between the transfer functions A;(/c) and powers of the power spectrum 
P{k), with cxi, Pi, ■Ji, Si corresponding to const., P{k), P{k)^/^, P{k)'^^^ respectively (refer to eqns (21) 
and (24).) Just as we did with the /nl estimator (80), we can sum the a^^'s from the Gaussian maps with 
these functions to create filtered maps, 

Ma{x,n) ^T.im^l{x)afJ^, M,{x,n) ^ J]7K^)aP^^, 

Im 

Mf,{x,n) ^EimPli^)<^?m^^ Ms{x,n) ^^dKx)ap^%P (114) 



Im 



Prom products of these maps in pixel space, we can now obtain explicit expressions for the nonGaussian 
a^'s in these two separable cases (compare with the bispectrum expressions (20) and (23)): 



_ local 



cquil 



J dxx'^ ^/3;(x) J d^hYi*^{h)Mo,ix,h) M0{x,h) + ^ai{x) J d^hYi*^{h)Mp{x,h) Mp{x,h 

2 J dxx^ -Wi{x) J d'^nYi*^{n) Maix,n) Mi3{x,n) - cxi{x) J d'^nYi'^{n} Mi3{x,n} Mp{x,h) 

-2diix) J d'^hYi*^{h)Msix,n)Msix,h) + 2^i{x) J d'^hYi*^{n) Mp{x,n) Msix,h) 

+2f3i{x) J d^nYil^{h)M^{x,h)Ms{x,h) + 2Si{x) J d^hYi*^{h)Mf^{x,i^^ . (115) 

In the top panel of fig. 12 we consider the contribution to the final C[^^ from the various terms ap- 
pearing in equation (115) taken separately. For example, we build a set of multipoles from the term 
J dxx'^^f3i{x) J d'^hYi^{h) Ma{x,h) Mi3{x,h) and compute the resulting power spectrum, neglecting all 
other terms, and so forth. We then compare the power spectra of the nonGaussian part to the input power 
spectrum of the Gaussian part for /nl = 100. Our procedure underlines what was pointed out in [32]: some 
terms in the separable approximations to both the local and equilateral shapes produce spurious divergences 
at low I's that are large enough to affect the final power spectrum of the map. More precisely, as can be 
seen in fig. 12, the biggest problems come from the terms f dxx'^'^Pi{x) f d^fi 1^*^.^(11) AIa{x, h) Aip{x, fi) and 
J dxx'^(3i(x) f d'^hYi'^{h) M^(x,h) Ms{x,h). In ref. [32], it was pointed out that the problem can be cir- 
cumvented for the local model by modifying the expression (115) of ajj^'^^ so as to eliminate the pathological 
term, while leaving the final bispectrum of the map preserved with a change of weight for the remaining 
term. The same approach can also be applied to the equilateral case, leaving new tailored expressions for 
the nonGaussian parts: 

= J dxx^ai{x) J d'^nYi^{n)M/3{x,n) M/3{x,n) , 

-3oci{x) [ d^hYC^{h)Mfi{x,h)Mp{x,n) -26i{x) [ d^hYi*^{h) Ms{x,n) Ms{x,n) 



„local 
^Im 



'^Im' = '^J dxx^ 



equil 



+27,(x) j d^uYil^{u)Mp{x,-h)Ms{x,-h) 



(116) 



It is easy to verify that these modified expressions produce the correct bispectra in the final maps, they 
are numerically stable and so allow the simulation of non-Gaussian maps of the local and equilateral type 
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Figure 12: Convergence properties of the standard separable functions used in the literature to represent local and equilateral 
models (top panel). Here the functions have been convolved with the transfer functions A; required in the bispectrum estimator 
or map-making algorithms. Note the poor scaling and divergence at low I for two of the separable combinations with the resulting 
power spectrum from non-Gaussianity rising to compete with the CMB power spectrum Ci's (/nl = 1). This poor scaling is 
contrasted with results for the tetrahedral polynomials qn{k) (lower panel). These remain bounded and roughly scale-invariant 
over the full multipole range, even for very high order polynomials. 



with given power spectrum and bispectrum. However, one can see how the necessity of looking at all 
the individual terms in the equations defining afj^, and the need to produce suitable modifications of the 
original formulae, means that the algorithm loses its generality. If additional shapes are considered then, in 
principle, different separation schemes could well encounter the problems outlined above. The good news is 
that the full generality of this approach is regained when the separation of the original shape is done using 
the eigenmode expansion introduced in this paper. 
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B. Map-making from arbitrary primordial shape functions 



In order to see why this happens, it is useful to write down the equation for aj^f in terms of our polynomial 
expression. Since we can decompose the three-point functions both at early and late times it is actually 
possible to generate a map in two different way. The closest method to the "standard" one, just outlined 
above, is the one that start from the early time decomposition. In this case the primordial shape 5(^1, k2, ks) 
is written as: 

Siki,k2,k3) = ^q;=Q„ = ^a^^^gp(fe2)gg(A;i)gr(fci) , (117) 

n pqr 

where the 271(^15^21^3) are formed from products of the tetrahedral polynomials qp{k) given in (54) and 
the a® are the coefficients of the eigenmode expansion for a given shape (recall the convenience of 

ordering the pqr with a single label n). The reduced angular bispectrum bij^i^ is obtained, as was shown 
in (76), by linearly projecting the primordial shape on the sphere using radiation transfer functions: 



^IInl^o:^ J x^dx 0y dki qp{ki) Ai^{ki) ji^{kix)^ (^/ "^^2 gr(fe2) Ai^(A;2) ji2(A:2X 

X j dks qsiks) Ai^{k3) ji^iksx)^ + 5 permutations 
AI/nl J^'^nJ x^dx qll,ql'qil , where g^(x) = dk qp(k) Ai{k) ji{kx) , 



(118) 



Substituting equation (118) into 110, and using the standard technique of decomposing the integrals into 
tractable products of one-dimensional integrals, after some algebra, we obtain the general expression for 



-.NG. 
hm ■ 



NG 
aim -^^ 



^ J dxx\'p{x) J d^hYr{n)M^{h,x)M^{h,x), (119) 



pqr*-*n 



where the M^(n, x) arc filtered maps found by summing a set of Gaussian a^'s with the convolved tetra- 
hedral polynomial (refer to eqn (80) 



M^{h,x)=^ql "Sr^ = T.lj Qp{k)Ai{k)ji{kx)dk 
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(120) 



The Mp{ii,x) and qp{x) are now the analogues of the Ma{x,n), Mg{x,h), Mj{x,n), Ms{x,n) and 
ai, /Sij'fi, di defined above. Note that it is not strictly necessary here to include cyclic permutations in 
(119) running over the indices {p,q,r}, as these are incorporated automatically. Further efficiencies can 
be achieved by exploiting the freedom to reorder terms in(119), taking out the polynomial qp of highest 
order and convolving the maps with the two lower order polynomials; this is not necessitated by stability 
requirements (see below). 

In principle, the numerical instabilities which cause problems for the standard separable approximations, 
could now affect the angular integrals over the polynomials qp given in (119). However, as shown in fig. 12, 
this is not the case. The key point is that all the functions qp{x) now scale as jij^ppj (see fig. 12), that is, 
in the same way as the non-pathological (3i{x) term in the standard local and equilateral decompositions. 
For this reason the spherical harmonic projection of a product of two Mp{h,x) maps is expected to have 
similar scaling properties as the term J d'^nYI^*{n)Mp{x, fi)^. This last integral was previously shown to be 
stable at low multipoles, as discussed for the local case in ref. [32]). Thus all the integrals in equation (119) 
are going to be well-behaved at low I's. Since the shape-dependent information is in the coefficients of the 
expansion a® and not in the precomputed qf{x) modes, we are able to produce numerically stable results 
for any possible shape. Numerical tests were carried out for the local and equilateral case, confirming the 
previous statements. We suggest, therefore, that the eigenmode expansion provides a numerically stable 
and efficient means by which to generalize the algorithm in ref. [31] to non-separable bispectrum shapes. 
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Figure 13: Simulated nonGaussian CMB maps Irom the equilateral model, created using the primordial map-making method 
(119). The upper panel shows a map simulation with /nl = 400 (barely discernable from the underlying Gaussian template), 
the middle panel shows a map with a large NG signal with /nl = 4000, while the lower panel shows the /nl = 400 case above 
in a WMAP-realistic context using the KQ75 mask and with inhomogeneous noise added. 
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C. Simulated maps from general CMB bispectra 

It is useful to recap the discussion above by using separable mode expansions to create simulated maps 
at late times from a given CMB power spectrum Ci and reduced bispectrum bij^^. As before with the /nl 
estimator, removal of the convolution with transfer functions, makes the late-time method much simpler 
and more transparent. We begin with the same expression 

"im = afm + fNLd^m ' (121) 

where 



O-lm 1^ Oi^kh^kmiW — — . (122j 

Ii,l2,mi,m2 ^ ^ 



Now we expand the CMB bispectra using the eigenmode decomposition using weight functions motivated 
by the estimator (refer to (89) 

^=^^.M3=E«"2- (123) 

where vi is a separable weight factor chosen to remove scaling from the CMB bispectrum, improving 
decomposition convergence. These weight factors are important for this late-time map-making method 
because they help remove the scaling of the ^/Cl term in the M^(n) filtered maps, making there power 
spectrum flatter (the analogue of the problem discussed above for primordial map simulations). We can 
rewrite the non-Gaussian part as 



= E^M} / dhYUn)M^{n)M^y{n) , (124) 

pqr<-^n 

where the M^(n) are defined in (92) and summed with Gaussian a^^s. 

This method is straightforward to implement for a given theoretical bij^i^ and it is highly efficent. For 
example, it can produce simulated maps in 64 seconds for I = 500 with 16 eigenmodes (6 polynomials). It 
has the advantage that, as it depends only on the CMB bispectra, it can also be used to simulate maps for 
bispectra produced by late time effects, like cosmic strings, gravitational lensing and secondary anisotropies. 
Plots of the non-Gaussian part of simulated maps can be seen in fig. 14 for the non-separable DBI inflation 
and cosmic string models. 



VI. DIRECT COMPARISON OF BISPECTRUM ESTIMATORS 



We have developed two complete numerical pipelines, implementing the eigenmode decomposition meth- 
ods described in the previous sections. For a generic primordial shape S{ki, k2, k^) or a given CMB bispec- 
trum bij^i^, an expansion in monomial symmetric polynomials Q„ is performed followed by the generation 
of nonGaussian map simulations. Bispectrum estimators are then appHed to the map simulations in order 
to verify that the input /nl can be properly recovered together with the expected variance. Both 'early 
time decomposition' and 'late time decomposition' /nl estimators have been fully implemented. The former 
starts from an expansion of the primordial shape S{ki, k2, k^) in Fourier space while the latter starts from 
an expansion of the reduced angular bispectrum bij^i^ in harmonic space, where in the second case radiation 
transfer functions have already been included in the expression for bij^i^. The redundancy provided by the 
two alternative pipelines provide a further check of the reliability of the final results. 

Since our purpose in this paper is to introduce the eigenmode expansion method, and to test its im- 
plementation, we will primarily apply our pipelines to map simulations, leaving detailed analysis of real 
datasets over a wider range of models for future publication [29] . As this is a proof-of-concept paper, we 
will mainly limit ourselves to the study of the simple equilateral family of models. This is because it is 
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Figure 14: Simulated maps for nonGaussian models using the late-time map-making method (124); this only includes the aj^„ 
contribution. The upper panel shows a non-Gaussian CMB map from cosmic strings obtained using the analytic expression 
for the string bispectrum (29). The lower panel shows a simulated nonGaussian map for an equilateral model. When added to 
its Gaussian counterpart map from af^ at an amplitude /nl = 600, this equilateral map was used for the bispectrum recovery 
illustrated in figs. 16, ??. Note the red colour cast from negative /nl and blue from positive. 

already well-studied in the literature (see e.g. [3, 25, 31]), which enables a useful comparison between the 
outcome of our numerical pipelines and previously published results for the equilateral shape. Moreover, 
the equilateral case does not require sophisticated noise analysis, unlike the local model. However, we will 
briefly consider other non-separable models outlined earlier in the introduction, such as the related DBI 
model and the cosmic string bispectrum. We note that from the point of view of the eigenmode decomposi- 
tion, the formal separability of the equilateral shape is irrelevant; it does not cause early termination of the 
expansion series which is nearly identical to the non-separable DBI model (see fig. 10). Having established 
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the reliability of the eigenmode expansion method here, in a forthcoming publication [29] we will apply it 
to the study of families of non-separable shapes using WMAP5 data. 

A. Simulated observational maps 

Using the algorithm described in section V, we generated a set of 100 equilateral CMB maps with both the 
primordial and late-time decomposition pipelines. We worked at roughly WMAP resolution with l^ax = 500 
and HEALpix nside = 512, corresponding to a pixel number Npix ~ 10^. We then applied both our 
primordial and late-time estimators to both our primordial and late-time sets of simulated maps in all 
combinations. We found that in all cases the map-making methods gave consistent results, producing 
simulated maps from which the correct /nl could be reliably recovered with the correct variance. Results 
for both primordial and late-time estimators on the same set of 50 equilateral maps (with and without the 
mask and inhomogeneous noise) with /jvl = 300 can be seen in fig. 15. We observe that the two estimators 
produce consistent results on the same maps. Of course, there is some small variation between the results as 
the two estimators can be regarded to be independent but this proved always to be well within the variance. 

In addition, we extracted the equilateral configurations Bm of the bispectrum from the maps and com- 
pared the average over all the simulations to the semi-analytic expectations obtained from the standard 
decomposition of the equilateral shape in terms of a;(x), Pi{x), 7;(x), di{x) (refer to eqns (21) and (24)). 
The recovered equilateral bispectrum values were in very good agreement between the semi-analytic predic- 
tion from the "standard" a,j3,j,S decomposition and the simulations, based on our eigenmode expansion, 
thus showing consistency with previous approaches. 

Finally, we reiterate that this general approach to map simulation was highly efficient, producing Planck 
resolution maps for the equilateral model on short timescales. This made estimator validation through 
Monte Carlo simulatoins easily achievable with only modest resources. For other well-behaved bispectra, 
such as the cosmic string model, the general method proved robust. Examples of non-separable maps 
already have been discussed and shown in fig. 14. 

B. Primordial and late-time /nl estimators 

Choosing an input value /nl = 300 for the sets of equilateral map simulations described above, we 
compared results from both the primordial and late-time bispectrum estimators. In order to verify the 
consistency of the two methods we selected the late-time map sets and applied both estimators to it. The 
tests were performed starting from a noiseless full-sky map and then more realistic simulations were used, 
including partial sky-coverage and an anisotropic noise component. The rms noise was obtained by coadding 
WMAP V and W channel using the same scheme as the one adopted for nonGaussian analysis by the WMAP 
team [3]. The sky-coverage was done using the KQ75 mask, also adopted by the WMAP5 team for their 
/nl analysis. Only the approximate form (??) of the estimator is used, and not the full form (33) including 
the full covariance matrix and a linear term. Note however that this approximation has been demonstrated 





Ideal simulations 


WMAP5 simulations 




Average 


St. Dev. 


Average 


St. Dev. 


Primordial estimator 


292.9 


104.8 


297.7 


152.1 


Late-time estimator 


300.6 


104.9 


278.7 


160 


Internal st. dev. 


38.5 


102.6 



Tabic I: Results obtained from the application of the primordial and latc-timc estimators as described in the text. In the first 
two columns, labeled by 'Ideal simulations', we consider ideal full-sky noiseless measurements, while in the last two columns, 
labeled by 'WMAP5 simulations' we include noise and sky coverage in order to simulate a WMAP5-realistic experiment (sec 
text for further explanation). We apply both estimators to a single set of maps, in this case created using the late-time mode 
expansion approach. In the last row, we calculate the difference between the /nl recovered by the two techniques, map by map 
for 100 maps, and report the final internal standard deviation between the methods. 
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Figure 15: Recovery of /nl from 50 simulated maps of the equilateral model, showing a direct map- by-map comparison between 
the primordial estimator (83) (blue) and the CMB estimator (95) (red). Ideal map recovery is shown in the top panel, while 
recovery for WMAP-realistic maps is shown below with beam, inhomogeneous noise and mask included (BNM). Both methods 
recovered the input /nl = 300 with a variance of approximately A/nl = 105 (clean) and 150 (BNM). Note the overall 
consistency of the two independent estimators with a significantly lower variance evident between the methods A/nl = 30 
(clean) and A/nl = 103 (BNM). 



in several previous studies to work well for equilateral shapes. Moreover, for our purposes the approximate 
nearly-optimal estimator is all we need since it contains all the dependence on the theoretical ansatz and 
thus all the dependence on our eigenmode expansion, which is the primary concern for this initial validation 
process. 

We compared the /tvl recovered from each map using the two methods, as well as the final averages 
and variances. The variances were compared to expectations from Fisher matrix forecasts obtained both 
from our eigenmode expansion and from the 'standard' cxi, Pi, 7^, Si decompostion of the equilateral 
shape used to date in other nonGaussian analysis. In all cases the results were internally consistent and in 
agreement with Fisher matrix expectations, as summarized in table (I). This led us to conclude that the 
eignemode expansion method appears to be a reliable way to produce non-Gaussian CMB simulations and 
/nl estimators for primordial models, whether separable or otherwise. 

Having verified the two estimator's performance on simulated equilateral maps we then applied both of 
them to the WMAP5 data, coadding the V and W channels as discussed above. The primordial estimator 
obtained the result —174 < /^x"^^ < 434, which is consistent with the existing constraints obtained using 
standard separable primordial approach (given the caveat that a number of these results have now been 
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Figure 16: Recovered spectral coefficients from the late time estimator (89) from a single map simulation for an equilateral 
model with /nl = 600 (or normalised relative to the local model _/<nl ~ HO); NG map simulation shown in the lower panel 
of fig. 14. In both panels, the original d;^ decomposition coefficients for the theoretical model are shown for comparison 
(blue). In the upper panel, the coefficients recovered from the single realisation are shown, with error bars {2a) estimated 
from 100 Gaussian maps. In the lower panel, the (3^ are recovered in a WMAP-realistic context using the KQ75 mask with 
inhomogeneous noise added. Note that the provide a remarkably good fit to thea^ given the significance of the non-Gaussian 
signal. 



superseded [10]). The constraint using the late-time estimator was -90 < /^'["'' < 550 which is shghtly 
larger than the primordial result but still well within the la range. As the late-time estimator can be 
regarded as independent of the primordial estimator, some variation is to be expected as we have seen 
already in fig. 15. The difference between the estimators is consistent with the internal variance of 103 
noted in table 1 for equilateral map simulations in a WMAP-realistic context. We conclude that both the 
primordial and late-time estimators appear to be performing up to expectation. 

Using the method described in eqn (101), we can endeavour to recover the full bispectrum from a given 
map. To illustrate this capability, we created a single map realization from an equilateral bispectra with 
/nl = 600, that is, a map with a 4a nonGaussian signal. We then used the late-time estimator to recover 
the and mode coefficients described in (89). Recall that for results of sufficient significance, the 
should approximate the original theoretical model coefficients a^, that is, those used to generate the 
simulated map. We estimated the variance in each of the eigenmodes by applying the same method to 100 
Gaussian simulations. The results for the orthornomal coefficients are plotted in fig. 16 for both ideal 
maps and for maps with inhomogeneous noise added and a mask applied. We see that we recover the first 7 
modes well from the ideal map but the results from the map containing noise and mask are somewhat less 
encouraging. Clearly, more work is required to control noise and mask effects at higher mode numbers. By 
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Figure 17: Recovered 3D bispectrum using the late-time mode decomposition method (101) from a single map simulation for 
an equilateral model with /nl = 600 (or relative to the local model, Jnl ~ HO); this figure shows the reconstruction of the 
bispectrum from the expansion modes illustrated in fig. 16. The left panel represents the original theoretical bispectrum 
used to construct a single realisation of the map (just like those shown in fig. 14). The middle panel represents the recovered 
bispectrum from the ideal map, while the right panel represents the recovery in a WMAP-realistic context using the KQ75 
mask with inhomogeneous noise added. The main feature of the bispectrum, that is, the primary acoustic peak appears to be 
evident even in the noisy cut-sky case, given the significant nonGaussian signal. 

plotting the 3D bispectra from the reconstructions, see figure 17, we observe that it is possible to recover 
the main accoustic peak and some basic features of the CMB bispectra. We will address the challenging 
issues associated with bispectrum reconstruction in greater detail elsewhere [29]. 

VII. CONCLUSIONS 

We have now implemented two comprehensive and independent pipelines for the analysis and estimation 
of general primordial or CMB bispectra. Both methods are based on dual mode expansions, exploiting a 
complete orthonormal eigenmode basis to efficiently decompose arbitrary bispectra into a separable poly- 
nomial expansion. These separable mode expansions, whether at late or early times, allow a reduction of 
the computational overhead to easily tractable levels, whether calculating the reduced bispectrum bi^i^i^, 
generating Planck resolution nonGaussian map simulations, or directly estimating /nl from simulations or 
real data sets. The method exploits the smoothness of the pattern of acoustic peaks observed in calcula- 
tions reviewing all well-behaved primordial models, implying the rapid convergence of the corresponding 
mode expansions. While most calculations are performed in the separable basis, a final rotation of mode 
coefficients into the orthonormal frame allows for a simple interpretation of the contributions to /nl using 
Parseval's theorem. In fact, the completeness of the orthonormal eigenmodes means, in principle, that it is 
straightforward to extract and reconstruct the full CMB bispectrum from the data, assuming the presence 
of a sufficiently significant nonGaussian signal. 

The main purpose of this paper has been to present a detailed theoretical framework for /nl estimation 
using separable eigenmode expansions, irrespective of the specific polynomials or other basis functions em- 
ployed. However, we have also presented some numerical results from the pipelines we have implemented, 
chiefly for the equilateral model where there are extensive published results for direct comparison. An im- 
portant milestone for the validation of this approach has been the development of a robust and reliable mode 
expansion method for generating map simulations from arbitrary bispectra. While generalising previous 
methods applied to specific separable cases, we noted that the scale-invariance of the polynomial expansion 
modes eliminates numerical instabilities that previously had to be circumvented on a case-by-case basis. 
Given convergent mode expansions for well-behaved bispectra, high resolution map simulations for a wide 
variety of models can easily and efficiently be generated, with several examples illustrated here including 
late-time cosmic strings. The many map simulations created for the equilateral model with both primordial 
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and late-time methods showed consistency in expected variance and /nl recovery. 

The primordial and late-time /nl estimators using mode expansions were tested successfully on the 
simulated equilateral maps, matching expectations for semi-analytic Fisher matrix forecasts and providing 
consistent unbiased results for /nl- This was achieved for both ideal maps and in a WMAP-realistic 
context, incorporating beams, anisotropic noise and a mask. Application of the estimators to the WMAP5 
data gave constraints on the equilateral model consistent with each other and previously published results. 
These encouraging results suggest that the approach will provide a robust and general framework for /nl 
estimation for the wide variety of non-separable models which remain to be constrained [11]. For single 
equilateral map simulation with /nl = 600, we were able to demonstrate a reasonable correspondence 
between the theoretical and recovered mode expansion coefficients, while also being able to recover key 
features of the full CMB bispectrum. However, a detailed discussion of such prospects has been left for 
a future publication [29]. We have also left aside for discussion elsewhere a more sophisticated treatment 
of sky cuts and inhomogeneous noise, which is more important for the analysis of the local model, as 
well as the potential for incorporating polarisation data. Challenges remain for the full implementation 
of the primordial and late-time pipelines at Planck resolution, but the generality and robustness of this 
methodology suggests that it should prove to be a useful tool for exploring and constraining a much wider 
class of nonGaussian models. 

VIII. ACKNOWLEDGEMENTS 

We are very grateful for informative discussions with Xingang Chen and Kendrick Smith and wc have 
also benefitted from useful conversations with Martin Bucher, Anthony Challinor, Olivier Forni, Enrique 
Martinez-Gonzalez and Bartjan van Tent. The ongoing development of these methods has been regularly 
reported at Planck Working Group 4 (NonGaussianity), where we have been grateful for feedback and to 
learn of related work. Simulations were performed on the COSMOS supercomputer (an Altix 4700) which 
is funded by STFC, HEFCE and SGI. JRF, ML and EPS were supported by STFC grant ST/F002998/1 
and the Centre for Theoretical Cosmology. 



[1] Juan Martin Maldacena. Non-Gaussian features of primordial fluctuations in single fleld inflationary models. JHEP, 
05:013, 2003. 

[2] Viviana Acquaviva, Nicola Bartolo, Sabino Matarrese, and Antonio Riotto. Second-order cosmological perturbations from 

inflation. Nud. Phys., B667:119 148, 2003. 
[3] E. Komatsu ct al. Five- Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Interpretation. 

Astrophys. J. SuppL, 180:330 376, 2009. 
[4] Amit P. S. Yadav and Benjamin D. Wandelt. Evidence of Primordial Non-Gaussianity (/nl) in the Wilkinson Microwave 

Anisotropy Probe 3- Year Data at 2.8a. Phys. Rev. Lett, 100:181301, 2008. 
[5] A. Curto et al. WMAP 5-year constraints on fnl with wavelets. 2008. 

[6] Oystein Rudjord et al. An Estimate of the Primordial Non-Gaussianity Parameter /nl Using the Needlet Bispectrum 

from WMAP. Astrophys. J., 701:369-376, 2009. 
[7] Kendrick M. Smith, Leonardo Senatore, and Matias Zaldarriaga. Optimal limits on /]vl^' from WMAP 5-year data. JCAP, 

0909:006, 2009. 

[8] ESA Planck Consortium. Planck: The Scientific Programme (Planck Blue Book). 2005. 

[9] Eiichiro Komatsu, David N. Spergel, and Benjamin D. Wandelt. Measuring primordial non-Gaussianity in the cosmic 

microwave background. Astrophys. J., 634:14-19, 2005. 
[10] Leonardo Senatore, Kendrick M. Smith, and Matias Zaldarriaga. Non-Gaussianities in Single Field Inflation and their 

Optimal Limits from the WMAP 5-year Data. 2009. 
[11] J. R. Fergusson and E. P. S. Shellard. The shape of primordial non-Gaussianity and the CMB bispectrum. Phys. Rev., 

D80:043510, 2009. 

[12] Mark Hindmarsh, Christophe Ringeval, and Teruaki Suyama. The CMB temperature bispectrum induced by cosmic 
strings. Phys. Rev., D80:083501, 2009. 

[13] D. M. Regan and E. P. S. Shellard. Cosmic String Power Spectrum, Bispectrum and Trispectrum. 2009. 
[14] J. R. Fergusson and Edward P. S. Shellard. Primordial non-Gaussianity and the CMB bispectrum. Phys. Rev., D76:083523, 
2007. 



41 



[15] Emiliano Scfusatti, Michele Liguori, Amit P. S. Yadav, Mark G. Jackson, and Enrico Pajer. Constraining Running Non- 

Gaussianity. 2009. 

[16] A. Curto, E. Martinez-Gonzalez, and R. B. Barreiro. Improved constraints on primordial non-Gaussianity for the Wilkinson 

Microwave Anisotropy Probe 5-yr data. Astrophys. J., 706:399-403, 2009. 
[17] Martin Bucher, Bartjan Van Tent, and Carla Sofia Carvalho. Detecting Bispectral Acoustic Oscillations from Inflation 

Using a New Flexible Estimator. 2009. 
[18] Daniel Babich, Paolo Creminelli, and Matias Zaidarriaga. The shape of non-Gaussianities. JCAP, 0408:009, 2004. 
[19] G. I. Rigopoulos, E. P. S. Shellard, and B. J. W. van Tent. Large non-Gaussianity in multiple-field inflation. Phys. Rev., 

D73:083522, 2006. 

[20] Xingang Chen and Yi Wang. Large non-Gaussianities with Intermediate Shapes from Quasi-Single Field Inflation. 2009. 

[21] Paolo Creminelli. On non-gaussianities in single-field inflation. JCAP, 0310:003, 2003. 

[22] Mohsen Alishahiha, Eva Silverstein, and David Tong. DBI in the sky. Phys. Rev., D70: 123505, 2004. 

[23] Xingang Chen, Min-xin Huang, Shamit Kachru, and Gary Shiu. Observational signatures and non-Gaussianities of general 

single field inflation. JCAP, 0701:002, 2007. 
[24] Daniel Babich. Optimal Estimation of Non-Gaussianity. Phys. Rev., D72:043003, 2005. 

[25] Paolo Creminelli, Alberto Nicolis, Leonardo Senatore, Max Tegmark, and Matias Zaidarriaga. Limits on non-Gaussianities 

from WMAP data. JCAP, 0605:004, 2006. 
[26] Amit P. S. Yadav et al. Fast Estimator of Primordial Non-Gaussianity from Temperature and Polarization Anisotropies 

in the Cosmic Microwave Background II: Partial Sky Coverage and Inhomogeneous Noise. Astrophys. J., 678:578, 2008. 
[27] E. Komatsu et al. Five- Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations :Cosmological Interpretation. 

2008. 

[28] Huiyuan Li and Yuan Xu. Discrete Fourier analysis on a dodecahedron and a tetrahedron. 2008. 

[29] James R. Fergusson, Michele Liguori, and Edward P. S. Shellard. CMB constraints on non-Gaussianity, in preparation. 
2009. 

[30] Ian G Moss and Chun Xiong. Non-gaussianity in fluctuations from warm inflation. JCAP, 0704:007, 2007. 

[31] Kendrick M. Smith and Matias Zaidarriaga. Algorithms for bispectra: forecasting, optimal analysis, and simulation. ArXiv 

Astrophysics e-prints, 2006. 

[32] Duncan Hanson, Kendrick M. Smith, Anthony Challinor, and Michele Liguori. CMB lensing and primordial non- 
Gaussianity. Phys. Rev., D80:083004, 2009. 



